In this paper we study the initial boundary value problem for the system −div (I + mm T )∇p = s(x), mt − α 2 ∆m + |m| 2(γ−1) m = β 2 m • ∇p∇p in two space dimensions. This problem has been proposed as a continuum model for biological transportation networks. The mathematical challenge is due to the presence of cubic nonlinearities, also known as trilinear forms, in the system. We obtain a weak solution (m, p) with both |∇p| and |∇m| being bounded. The result immediately triggers a bootstrap argument which can yield higher regularity for the weak solution. This is achieved by deriving an equation for (I + mm T )∇p • ∇p j , j ≥ 1, and then suitably applying the De Giorge iteration method to the equation.