Global existence of strong solutions to a biological network formulation model in $2+1$ dimensions

In this paper we study the initial boundary value problem for the system $-\mbox{{div}}\left[(I+\mathbf{m} \mathbf{m}^T)\nabla p\right]=s(x),\ \ \mathbf{m}_t-\alpha^2\Delta\mathbf{m}+|\mathbf{m}|^{2(\gamma-1)}\mathbf{m}=\beta^2\mathbf{m}\cdot\nabla p\nabla 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 $(\mathbf{m},p) $ with both $|\nabla p|$ and $|\nabla\mathbf{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 $\left((I+\mathbf{m} \mathbf{m}^T)\nabla p\cdot\nabla p\right)^j, j\geq 1$, and then suitably applying the De Giorge iteration method to the equation.


Introduction
Continuum models for biological transportation networks have received tremendous attention recently. We refer the reader to [2] for a rather comprehensive survey of the subject. The most well known model is the one proposed by Hu and Cai [8,9]. It deals with the scalar pressure function p = p(x, t) and the vector-valued conductance function m = m(x, t). They satisfy the system −div (I + mm T )∇p = s(x) in Ω T ≡ Ω × (0, T ), (1.1) m t − α 2 ∆m + |m| 2(γ−1) m = β 2 m · ∇p∇p in Ω T , (1.2) coupled with the initial boundary conditions p = m = 0 on Σ T ≡ ∂Ω × (0, T ), (1.3) m(x, 0) = m 0 (x) on Ω, (1.4) where Ω is a bounded domain in R N with boundary ∂Ω and T > 0.
The proof in [6] was based upon the following a priori estimates 1 2 Ω |m(x, τ )| 2 dx + D where τ ∈ (0, T ], Ω τ = Ω × (0, τ ), and p 0 is the solution of the boundary value problem −div[(I + m 0 m T 0 )∇p 0 ] = s(x), in Ω, (1.8) p 0 = 0 on ∂Ω. (1.9) Finite time extinction or break-down of solutions in the spatially one-dimensional setting for certain ranges of the relaxation exponent γ was carefully studied in [7]. Further modeling analysis and numerical results can be found in [1]. We also mention that the question of existence in the case where γ = 1 2 is addressed in [7]. In this case the term |m| 2(γ−1) m is not continuous at m = 0. It must be replaced by the following function However, the general regularity theory remains fundamentally incomplete. In particular, it is not known whether or not weak solutions develop singularities in 2 or higher dimensions. When the space dimension N is three, the initial value problem for the system (1.1)-(1.2) has been studied in [11], where the local existence of a strong solution and global existence of such a solution for small data are established. In addition, the author obtained a condition which a strong solution must satisfy if it blew up in finite time. However, the author specifically mentioned that his method there was not applicable to the case where N = 1 or 2. If N = 2, the same initial value problem was considered in [15]. Here the authors obtained a similar blow-up criterion to that in [11] and the global existence of a strong solution under the additional assumptions that α is sufficiently large and γ ≥ 1. As for the initial-boundary value problem for (1.1) and (1.2), Jian-Guo Liu and the author [12] obtained a partial regularity theorem for (1.1)-(1.4). It states that the parabolic Hausdorff dimension of the set of singular points can not exceed N , provided that N ≤ 3. A different form of partial regularity is obtained in [18]. If N = 2, then it is shown in [17] that p is continuous in the space variables and (m, p) are classical if they are stationary.
Of course, this theorem can start a bootstrap argument which results in even higher regularity. In fact, our result here implies the existence of the so-called strong solution [11,18] in the following sense. Indeed, it immediately follows that β 2 ∇p · ∇p − |m| 2(γ−1) m ∈ L ∞ (Ω T ). The classical regularity theory for the heat equation ( [10], Chapter IV) asserts that |m t |, |∆m| 2 ∈ L s (Ω T ) for each s ≥ 1. Thus we can write (1.1) in the form (1.10) tr (I + mm T )∇ 2 p + div(I + mm T )∇p = −s(x).
This enables us to conclude from the classical regularity theory for elliptic equations that p ∈ L ∞ (0, T ; W 2,s ) for each s ≥ 1. Moreover, we can differentiate this equation with respect x i to derive an equation for p x i . This combined with our high regularity assumption on the boundary ∂Ω ensures that we have p ∈ W 3,s (Ω) for some s > 1. Note that we only assume that γ > 1 2 . As a result, we cannot differentiate the term |m| 2(γ−1) m. This will prevent us from obtaining any estimates for the third order partial derivatives of m.
We follow [6] for the notion of a weak solution.
Note that the elliptic coefficients in (1.1) satisfy That is, (1.1) is only singular. This enables us to show that p is bounded [12]. In fact, we have p ∈ L ∞ (0, T ; C loc (Ω)) [17]. Unfortunately, this is not enough to trigger a bootstrap argument. We must have the Hölder continuity of p in the space variables to obtain the boundedness of m (see Lemma 2.2 below). Instead of trying to bridge this gap, we directly go after the boundedness of ∇p. This is motivated by a result in [3] where the author considered an elliptic equation of the form Here we have employed the Einstein summation. That is, repeated indices are implicitly summed over. An equation for ln (A∇u · ∇u) was derived to study critical points of u. In our case if we let we can derive an equation for ψ ≡ (A∇p · ∇p) j for each j ≥ 1. To be specific, we have for some functions H, h, K. In particular, H, h, K are only bounded by |m|, |∇m|. The trade-off is that equation (1.13) is both degenerate and singular. We overcome these singularities by suitably modifying the classical De Giorge iteration method. Even though the derivation of (1.13) is inspired by a result in [3], there are some major differences. The most prominent one is that we have not been able to impose the normalization condition a 11 a 22 − a 2 12 = 1 as did in [3]. Doing so would have changed the smallest eigenvalue of the coefficient matrix to 1 √ 1+|m| 2 , which is not bounded away from 0 below because we do not have the a priori knowledge that m is bounded. The resulting estimate for ψ would be useless to us. As we shall see, not being able to normalize the coefficient matrix causes many complications. Also, if we just have an equation for ln (A∇p · ∇p) as did in [3] or p x i , i = 1, 2, that will not help us at all. The fact that j in (1.13) can be arbitrarily large is crucial to the success of our argument.
We do have local existence of strong solutions [18]. Our earlier discussion suggests that such solutions satisfy the above condition. Thus our result simply says that strong solutions do not blow up in finite time. We also remark that our assumption on s(x) is more than necessary. If we check our calculations carefully, it is enough to assume s(x) ∈ L q (Ω) for q sufficiently large. Then we must use functions in L ∞ (Ω) ∩ W 1,q (Ω) to approximate s(x) in order to gain enough regularity for p to justify the calculations in the derivation of (1.13). The rest of the paper is organized as follows. Section 2 is largely preparatory. Here we collect some relevant known results for later use. To justify all the calculations in Section 3, p must be sufficiently regular. At the end of this section we offer an approximation scheme which produce sufficiently regular solutions. In Section 3 we first derive (1.13). Then the proof of the main theorem is achieved in two stages. First we show that ∇p ∞,Ω T is bounded by ∇m ∞,Ω T . This is done via the De Giorge iteration method. Then we prove that ∇m ∞,Ω T is also bounded ∇p ∞,Ω T . The key to the success of our argument is that j in the definition of ψ can be arbitrarily large.
Let us make some remarks about notations. The letter c is used to denote a positive number whose value can be computed from given data. The capital letters such as A, B, · · · are often used to represent 2 × 2 matrices. The ij-entry of A is denoted by a ij . The boldface letters are used to denote vector quantities. The i-th entry of F is f i .

Preliminary results
In this section we first collect some formulas about differentiating matrix-valued functions. Then we prove that local Hölder continuity of p in the space variables implies the local boundedness of m. At the end of the section we present an approximation scheme which gives necessary regularity to justify the calculations in the subsequent section.
Denote by M 2×2 the space of all 2 × 2 matrices. We invoke the following notation conventions Denote by ∇ 2 p the Hessian of p. Then we have The following identities will be frequently used We also need the interpolation inequality The next lemma deals with sequences of nonnegative numbers which satisfy certain recursive inequalities.
In fact, by Theorem 7.15 in ( [5], p. 162), there is a positive number c 0 such that Fix y ∈ Ω. For r ∈ (0, dist(y, ∂Ω)) we choose a smooth cutoff function ξ with the properties We use (p − p y,r (t))ξ 2 as a test function in (1.1) to get Br(y) By choosing s sufficiently close to 1, we can find a positive number ε such that Take the dot product of (1.2) with m to obtain Consider the problem By the comparison principle, we have The right-hand side term in (2.22) satisfies (2.19), a result in [19] asserts that w is Hölder continuous. This implies the desired result. The proof is complete.

Boundedness for ∇p and ∇m
In this section we will offer the proof of the main theorem. We shall begin by deriving (1.13).
Let A be given as in (1.12). Recall from (2.7) that We can write (1.1) in the form As in [3], we introduce the following functions v = A∇p · ∇p, div(A∇φ) = H · ∇φ + h + divK.
for some functions H, h, K.
We will identify the three functions in the above equation at the end of the proof.
Proof. The identity in [3] is still valid here. To see this, we compute from (2.6) and (2.8) that The last step is due to the fact that A is symmetric. The first two terms on the right-hand side of the above equation are troubling. One contains third order partial derivatives of p, while the other is quadratic in E. It turns out that both terms can be represented in terms of det(∇ 2 p). After we substitute them back into (3.11), the det(∇ 2 p) terms get canceled out. We shall do this by finding a suitable formula for the matrix D defined by An elementary calculation shows that the four entries of D are as follows Using (3.2), we obtain Now we are in a position to calculate that div(AE) = 2div(A∇ 2 pA∇p) = 2∇ 2 p : (A∇ 2 pA) + 2div(A∇ 2 pA)∇p. (3.17) Applying the formula for A∇ 2 pA yields Here we have used (2.9) and the fact that div Collecting the preceding two results in (3.17) gives As for AE · E, we have We are ready to calculate We still need to eliminate the second partial derivatives of p in the last term of of the preceding equation. If det(A) had been 1, then this term would be zero, and hence the proof would conclude. Since we do not have the benefit, we need to continue. We deduce from (3.9) and (3.2) that 2(a 11 p x 1 + a 12 p x 2 )p x 1 x 1 + 2(a 21 p x 1 + a 22 p x 2 )p x 1 x 2 = e 1 , (3.23) 2(a 11 p x 1 + a 12 p x 2 )p x 1 x 2 + 2(a 21 p x 1 + a 22 p x 2 )p x 2 x 2 = e 2 , (3.24) a 11 p x 1 x 1 + 2a 12 p x 1 x 2 + a 22 p x 2 x 2 = w. (3.25) Denote by E the coefficient matrix of the above system. Then By Cramer's rule, we have [−a 11 (a 11 p x 1 + a 12 p x 2 )e 1 − a 12 (a 11 p x 1 + a 12 p x 2 )e 2 ] This yields where A 1 = a 11 (a 11 p x 1 + a 12 p x 2 ) a 12 (a 11 p x 1 + a 12 p x 2 ) a 11 (a 12 p x 1 + a 22 p x 2 ) a 22 (a 11 p x 1 + a 12 p x 2 ) , A 2 = a 11 (a 12 p x 1 + a 22 p x 2 ) a 22 (a 11 p x 1 + a 12 p x 2 ) −(a 22 a 11 − 2a 2 12 )p x 1 + a 12 a 22 p x 2 a 22 (a 12 p x 1 + a 22 p x 2 ) , A 3 = −(a 11 p x 1 + a 12 p x 2 ) 2 −(a 11 p x 1 + a 12 p x 2 )(a 12 p x 1 + a 22 p x 2 ) −(a 11 p x 1 + a 12 p x 2 )(a 12 p x 1 + a 22 p x 2 ) −(a 12 p x 1 + a 22 p x 2 ) 2 .
In summary, we have The proof is finished by setting We would like to remark that the last part in our proof only works for two space dimensions. If the space dimension had been three, we would have six second order partial derivatives. But (3.9) and (3.2) would only give us four equations. Thus the same argument would fail. However, in the context of our proof, the last part becomes necessary only because we cannot normalize the coefficient matrix. Even if we could have done this, our argument would still only work for the two space dimensions. As we can easily see, it is not possible to represent AE · E in terms of det ∇ 2 p if the space dimensions are bigger than or equal three.
It immediately follows that With these in mind, we can derive that |G| ≤ c|m||∇m|, (3.37) |w| ≤ c|m||∇m||∇p| + |s(x)|, (3.38) We can easily deduce that Let j be given as in the theorem and define The equation satisfied by ψ is given by (1.13). Now fix a point x 0 ∈ Ω. Then pick a number R from (0, dist(x 0 , ∂Ω)). Define a sequence of concentric balls B Rn (x 0 ) in Ω as follows: Choose a sequence of smooth functions θ n so that for each x ∈ R 2 , and (3.53) Hence, (3.57) K n ≥ 1 for each n.
We use θ 2 n+1 (ψ − K n+1 ) + as a test function in (1.13) to obtain This together with (3.34) and (3.58) implies The last term in (3.61) can be estimated as follows: With this in mind, we estimate Here χ S n+1 (t) is the indicator function of the set S n+1 (t). Similarly, Plugging the preceding results into (3.61), we obtain We pick a number r from the interval (1, ∞). Define We conclude from (3.69) that . By Poincaré's inequality, we have We easily see that Substituting this into (3.73) yields In view of Lemma 2.1 and (3.55), it is enough for us to take Similarly, by (H1), (3.46), and (3.47), we have Collecting the preceding estimates in (3.77) and taking the j th root of the resulting inequality, we arrive at By an argument in ( [5], p. 303), we can extend the above estimate to the whole Ω. That is, we have We are ready to prove the main theorem.
It is important to note that for our argument to work we must be able to choose j big enough.