We consider n nonintersecting Brownian motion paths with p prescribed starting positions at time t=0 and q prescribed ending positions at time t=1. The positions of the paths at any intermediate time are a determinantal point process, which in the case p=1 is equivalent to the eigenvalue distribution of a random matrix from the Gaussian unitary ensemble with external source. For general p and q, we show that if a temperature parameter is sufficiently small, then the distribution of the Brownian paths is characterized in the large n limit by a vector equilibrium problem with an interaction matrix that is based on a bipartite planar graph. Our proof is based on a steepest descent analysis of an associated (p+q)×(p+q) matrix-valued Riemann–Hilbert problem whose solution is built out of multiple orthogonal polynomials. A new feature of the steepest descent analysis is a systematic opening of a large number of global lenses.