Mathematical Analysis of a Kinetic Model for Cell Movement in Network Tissues

Mesenchymal motion describes the movement of cells in biological tissues formed by fiber networks. An important example is the migration of tumor cells through collagen networks during the process of metastasis formation. We investigate the mesenchymal motion model proposed by T. Hillen (J. Math. Biol. 53:585-616, 2006) in higher dimensions. We formulate the problem as an evolution equation in a Banach space of measure-valued functions and use methods from semigroup theory to show the global existence of classical solutions. We investigate steady states of the model and show that patterns of network type exist as steady states. For the case of constant fiber distribution, we find an explicit solution and we prove the convergence to the parabolic limit.


(Communicated by Kevin Painter)
Abstract. Mesenchymal motion describes the movement of cells in biological tissues formed by fibre networks. An important example is the migration of tumour cells through collagen networks during the process of metastasis formation. We investigate the mesenchymal motion model proposed by T. Hillen in [14] in higher dimensions. We formulate the problem as an evolution equation in a Banach space of measure-valued functions and use methods from semigroup theory to show the global existence of mild and classical solutions. We investigate steady states of the model and show that patterns of network type exist as steady states. For the case of constant fibre distribution, we find an explicit solution and we prove the convergence to the parabolic limit.
1. Introduction. Friedl and collaborators [11] observed mesenchymal tumour cells as they move in a field of collagen fibres and change their velocities according to the local orientation of the fibres. At the same time, the cells also remodel the fibres, primarily through expression of matrix-degrading enzymes (proteases) that cut selected fibres. In [14], the author introduced a mathematical model for this process of mesenchymal cell movement in fibrous tissues. Recent analysis of this and similar models [14,21,4,5,25] revealed the existence of biologically meaningful measure valued solutions, which correspond to tissue and cell alignment. Hence a sophisticated existence theory is needed. In this paper we will formulate the mesenchymal transport model proposed in [14] as a semilinear evolution equation in a Banach space of measure-valued functions. We apply classical theory of semigroups of operators and a Banach Fixed Point argument to show well-posedness of the problem (Section 3.1). With the correct theoretical framework in place, we are 1056 THOMAS HILLEN, PETER HINOW AND ZHI-AN WANG then able to classify possible steady states, whereby we introduce a new notation of pointwise steady states, which are meant to resemble the network patterns which were observed numerically in [21]. Moreover, we rigorously study the parabolic limit (diffusion limit) of the kinetic model in the measure-valued context. We show convergence to the diffusion limit for constant fibre distribution in Section 5.1.
The existence theory here employs a mild solution formulation which is based on a variation of constant formula. The solutions are functions in L 1 in space and measures in velocity. It turns out that this definition is too "weak" in the sense that it does not provide a nice representation of the global network patterns observed numerically. Hence here we introduce a sub-class of steady states, which we call pointwise steady states. First of all, we show that pointwise steady states do exist. Secondly, pointwise steady states allow for a representation of network patterns. Our results include a classification of possible network intersections.
In the model proposed in [14], undirected and directed tissues were distinguished. In undirected tissues (e.g. collagen), fibres are symmetric and both directions are identical, a situation that somewhat resembles a nematic liquid crystal [24]. In directed tissues, fibres are asymmetric and the two ends can be distinguished. From the mathematical point of view, which we adopt in the present paper, both cases are completely analogous. Hence we focus on the case of undirected tissues. We refer to [14] for the biological assumptions and the detailed mathematical derivation of the model.
The model studied here is specifically designed for mesenchymal cell movement in network tissues via contact guidance and degradation of the extracellular matrix (ECM). Painter [21] has extended this model in various directions. His model variations allow (i) to choose between amoeboid and mesenchymal motion, (ii) to place different weights between diffusive movement and movement by contact guidance, (iii) to include ECM degradation as well as production, (iv) to include ECM remodeling or lack thereof, (v) to study focussed protease release at the cell tip versus unfocussed ECM degradation via a diffusible proteolytic enzyme. Many of these modifications lead to the same pattern formation properties as observed for the initial model. All of these modifications show the same mathematical challenges, namely the description of aligned tissue as weak solutions and orientation driven instabilities. Hence we believe that the results which we present here are representative for a large class of kinetic models for cell movement in tissues and they can be generalized to many other cases.
In [14], the techniques of moment closure, parabolic and hydrodynamic scaling were used to study the macroscopic limits of the system that we later restate in equation (1). The resulting macroscopic models have the form of drift-diffusion equations where the mean drift velocity is given by the mean orientation of the tissue and the diffusion tensor is given by the variance-covariance matrix of the tissue orientations. Model (1) has been extended in [4,5] to include cell-cell interactions and chemotaxis. The corresponding diffusion limit was formally obtained in these papers.
In the case of chemotaxis, a system of a transport equation for the cell motion coupled to a parabolic or elliptic equation for the chemical signal was studied by Alt [1], Chalub et al. [3] and Hwang et al. [15,16]. Local and global existence of solutions were studied and the macroscopic limits were proved rigorously in [3,15,16]. However, these authors assumed the existence of an equilibrium velocity distribution for cells that is in L ∞ (V ) where V denotes the space of velocities. For the mesenchymal motion model here, it is necessary to allow for complete alignments of either fibres or cells, corresponding to Dirac measures on V or the space of directions, the unit sphere S n−1 . In particular, assumption (A0) in paper [3] does not apply here and hence their respective results can not be applied directly to the mesenchymal motion model.
In Section 2 we formulate the model and we introduce suitable function spaces and operators. Our first main result on global existence of measure-valued solutions is given in Section 3. In Section 4 we present a definition and classification of pointwise steady states. In Section 5 we assume that the fibre density q is a given function of x, t. In that case we find an explicit solution of the kinetic equation using the methods of characteristics. If moreover, the fibre distribution is constant in time and space, then we prove the convergence to a parabolic limit. It appears to be impossible to prove convergence to the parabolic limit for arbitrary time and space dependent fibre distributions. This confirms numerical observations of Painter [21], who investigated the mesenchymal motion model and found interesting cases of pattern formation of network type (see Figure 2). In the diffusion limit, however, the patterns disappear in the numerical simulation. This indicates that there is a significant difference in the asymptotics of the kinetic model and the diffusion limit for timely varying tissue networks.
2. Formulation of the problem.
2.1. The model. We briefly recall the kinetic model for mesenchymal motion from [14] for the undirected case. The distribution p(x, t, v) describes the cell density at time t ≥ 0, location x ∈ R n and velocity v ∈ V . Throughout the paper we assume that V is a product V = [s 1 , s 2 ] × S n−1 , where 0 ≤ s 1 ≤ s 2 < ∞ is the range of possible speeds. If s 1 = s 2 then we assume s 1 > 0. The fibre network is described by the distribution q(x, t, θ) with θ ∈ S n−1 , the (n − 1)-dimensional unit sphere in R n . A schematic of the model is given in Figure 1.  The model for mesenchymal motion from [14] reads where µ and κ are positive constants. The transport term v · ∇p indicates that cells move with their velocity. The right hand side of the first equation describes the reorientation of the cells in the field of fibres. Turning away from their old direction at rate µ, they turn into a new direction with a probability that corresponds to the fibre distribution q. The new speed is chosen from the interval [s 1 , s 2 ]. The cells degrade (at rate κ) those fibres that they meet at an approximately right angle while they leave fibres that are parallel to their own orientation unchanged. The exact definitions of the corresponding terms in system (1) requires some mathematical details that are given in the next subsection. The expressionsp,q, Π u (p) and A u (p, q) are defined in equations (2), (3), (5) and (6), respectively. Painter showed in [21] that the second equation of (1) arises if instead of ECM degradation one assumes that the cells realign the tissue. This would be the case for fibroblasts, who do remodel the fibre newtork without destroying it. In that case the term −A u measures the fibre degradation while Π u describes the fibre production such that the total amount of fibre mass is preserved.

Spaces and operators.
We show in Section 4 that Dirac measures occur as meaningful steady states. Hence we need to construct a solution framework that allows for measure-valued solutions. Let Ω = R n be the spatial domain in which particles are able to move.
Let B(V ) denote the space of regular signed real-valued (finite) Borel measures on V . For p ∈ B(V ) let p = p + − p − be its Hahn-Jordan decomposition and |p| = p + + p − its variation [6]. When equipped with the total variation norm (the following notations are used interchangeably throughout the paper) [6,Proposition 4.1.7]. Analogously, B(S n−1 ) will denote the Banach space of regular signed Borel measures on S n−1 equipped with the total variation norm. Naturally, we are interested in solutions taking values among non-negative measures only. Let We denote the positive cones of the spaces X 1 , X 2 and X by X + 1 , X + 2 and X + , respectively. We will write for those p ∈ X 1 for which the essential supremum is finite.
We define the following operators • The spatial mass density of a velocity distribution, Clearly, the operator¯is Lipschitz continuous. • The lifting of a measure on S n−1 to a measure on V , where m is a probability measure on [s 1 , s 2 ]. If s 1 = s 2 , then˜just maps a measure on S n−1 to the same measure on {s 1 }×S n−1 . In the paper [14] it was taken to be the normalized Lebesgue measure on [s 1 , s 2 ], which corresponds to the weight parameter ω defined in [14, equation (4)]. The choice m([s 1 , s 2 ]) = 1 guarantees that In particular, a function that takes values among the probability measures on B(S n−1 ) + is mapped to a function taking values among probability measures on B(V ) + . Since˜is a linear operator it is Lipschitz continuous. Additionally, we use the lifting to connect the measures on V and on S n−1 in a natural way as • The mean projection operator (for undirected fibres) For sake of simpler notation and to avoid difficulties whenp = 0, we introduce the operator Λ : Notice that Λ is linear and if ||p|| ∞ < ∞ then For sake of completeness we also state the directed version of the operator Λ, As said above, existence of solutions is shown completely analogously in the two cases. • The relative alignment operator again, using the notation from [14] A u (p, q) = Similarly to the introduction of Λ, we will work with Notice that B is bilinear and if ||p|| ∞ < ∞, then The operators Λ and B are Lipschitz continuous on bounded subsets. Let µ > 0 denote the turning rate and κ > 0 denote the rate of fibre degradation. The model (1) can be written as equality of measures 3. Existence results. To provide a framework for local and global existence of solutions we define Here where the integrals are Bochner integrals taking values in B(V ). Observe that the domain D(A) is dense in X, as it contains the space C ∞ (R n , B(V ))×X 2 of infinitely differentiable functions, which is dense in X [19, Theorem 2.16]. The operator A with domain D(A) is the generator of a positive C 0 -semigroup U(t) on the Banach space X (see also Theorem 1 in [2]).
Notice that the operator −v · ∇ is the collisionless transport operator occurring in the linear Boltzmann equation which has been studied by many authors, see [13,10], [18,Chapter 13] and the references therein. It generates a semigroup (in fact, a group) U 1 on the space X 1 via for Borel sets A ⊂ V . Clearly, the positive cone X + 1 is invariant under U 1 . The group U 1 preserves the L 1 -norm while for ||p 0 || ∞ < ∞ we have We denote the semigroup on X generated by the operator A from equation (8) by (U(t)) t≥0 . It has a diagonal structure where I denotes the identity on X 2 . In the operator norm, U satisfies ||U(t)|| L(X) ≤ 1 and for ||u 0 || ∞ < ∞ we obtain For a pair u = (p, q) ∈ X define the map || · || ∞ : X → [0, ∞] by and set D = {u ∈ X : ||u|| ∞ < ∞}. For every r > 0 the set D r = {u ∈ X : ||u|| ∞ ≤ r} is closed in X (in particular, the projection of D onto X 1 is closed with respect to the L 1 -norm on X 1 ). Indeed, let p n ∈ X 1 be a Cauchy sequence with ||p n || ∞ ≤ r for all n. Since X 1 is complete, it has a limit p. We claim that ||p|| ∞ ≤ r. Suppose that this were not the case, then there would be an ε > 0 and a set A ⊂ R n with Lebesgue measure |A| > 0 such that ||p(x)|| B(V ) ≥ r + ε for all x ∈ A. But then clearly the L 1 -norm would satisfy ||p n − p|| X1 ≥ ε|A| > 0, which is a contradiction.
) is continuous and u satisfies the integral equation where U(t) is the semigroup defined in equation (10). We call a function u = (p, q) : , and (ii) equation (12) holds.

Our first result is
Theorem 3.2. Assume that q 0 (x, S n−1 ) = 1 for almost every x ∈ R n , then the problem (12) has a unique global positive mild solution for every u 0 ∈ D ∩ X + .
3.1. Proof of Theorem 3.2. The proof of Theorem 3.2 is established in the following Lemmas.
Lemma 3.3. The right hand side of equation (7) defines a nonlinear map F : D → D, which maps D into itself The map F is Lipschitz continuous on bounded subsets of D.
Proof. Observe that for (p, q) ∈ D the productpq is well defined and in particular, For functions ϕ ∈ L ∞ (S n−1 ) and measures q ∈ B(S n−1 ) we define the product ϕq ∈ B(S n−1 ) by way of where M ⊂ S n−1 is a Borel set. This multiplication extends to functions in L ∞ (R n × S n−1 ) and L ∞ (R n , B(S n−1 )) and we have showing that F 2 takes values in X 2 . Computations similar to those just carried out give the local Lipschitz continuity of F on bounded subsets of D. For example, for We omit the remaining calculations.
Lemma 3.4. Equation (12) has a unique local mild solution that remains positive for u 0 ∈ X + .
Proof. We set up a Banach's Fixed Point argument, but we cannot work on D directly since that set is not complete. Hence we work with D R for some R large enough. For given u 0 ∈ D and fixed R, T > 0 we define This set E R,T is a complete metric space, with the metric given by For a function u ∈ E R,T we define this is again an element of C([0, T ], D) with Gu(0) = u 0 (since D is invariant under both the semigroup U and the nonlinearity F ). We have for u, v ∈ E R,T and by choosing T sufficiently small, it can be achieved that G is a contraction on the space E R,T . If v ∈ D R , then we have (see the proof of Lemma 3.3) Let u ∈ E R,T . We can estimate equation (15) By choosing R > 2||u 0 || ∞ (large) and T small, namely we can achieve that for all u ∈ E R,T . Hence the contraction G maps the complete metric space E R,T into itself and hence has a unique fixed point by the Banach Fixed Point Theorem. The positivity of solutions follows from the fact that the nonlinearity F is of multiplicative type. If either p or q becomes zero on a set at some time, the left hand side of equation (12) is non-negative.
Concerning the global existence of solutions, if T max (u 0 ) < ∞, then by [22] lim t Tmax(u0) However, in our system (7) a blow-up in finite time cannot occur as the following lemma shows.
Lemma 3.5. Let (p, q)(t) be a mild solution of equation (7) (equivalently, of (12)) taking values in D ∩ X + . Then for all t ∈ [0, T max ) and almost every x ∈ R n we have q(x, t, S n−1 ) = 1, and there exists a constant C > 0 such that Proof. Let (p, q) be a mild solution of equation (7) taking values in X + . The second component of equation (13) reads where I denotes the identity. We evaluate this relation at S n−1 and use the fact that We obtain We apply Gronwall's lemma and obtain The integrand is positive and bounded, hence by the assumption on q 0 we get q(x, t, S n−1 ) = 1 (16) for almost all x ∈ R n . For the L 1 -norm of p we notice first that since p is positive, it satisfies We evaluate the first equation of (7) and obtain as a consequence of (16) Integrating this equation over R n gives by the divergence theorem. For the L ∞ -part of p we use that F 1 (0, q) = 0 and the following fact We estimate from (13) This inequality warrants application of Gronwall's lemma 4. Steady states. In numerical simulations by Painter [21], shown in Figure 2, we find interesting network patterns which form from random initial data. Numerically, these patterns do not change after they have been established. We expect that the system (7) allows for these network patterns as steady states. In this section we will develop a theory of pointwise steady states which are candidates for the observed network patterns.
To describe steady states of (7) we introduce the bilinear turning operator Observe that in contrast to the paper of Chalub et al. [3] the turning kernel does not depend explicitly on v , i.e., the cells are reoriented regardless of their original orientation. For p ∈ ker L[q], we have p =qp.
Hence the orientation of the cells in a steady state is entirely given by the fibre distribution q. This reflects the fact that a perfect alignment of the cells with the underlying fibre network and only such a perfect alignment remains invariant under the turning operator L. The trivial steady state is a uniform distribution of fibres and cells: where the small bars indicate the mean direction and the gray color describes the degree of alignment. Light gray indicates highly aligned tissue, whereas dark gray/black indicates close to uniform distribution of directions. The simulations were done by K. Painter, and are described in detail in [21]. We are grateful to K. Painter who allowed us to use this figure for illustrative purposes.
is a steady state of (7) in L ∞ (R n , B(V ) × B(S n−1 )). The only steady state of this type in D ∩ X is obtained for = 0.
Proof. If q = dθ/|S n−1 | and p = q, thenp = and p =pq. The right hand side of the first equation of (7) is zero. For the second equation, we need to compute Λ(p) − B(p, q). We have for a β ≥ 0, which is independent of θ and x. We obtain Hence Λ(p) − B(p, q) = 0 and the right hand side of the second equation of (7) is zero as well. Notice that the measures dv and dθ are coupled in a natural way through (4).
To find other steady states, we need a weak formulation.
Definition 4.2. We say that (p, q) ∈ D ∩ X is a weak steady state of (7), if for each pair of test functions φ ∈ W 1,1 (R n , C(V )), ψ ∈ C 0 (R n , C(S n−1 )), Notice that in this definition, φ and ψ are real-valued functions in the variables v and θ, respectively, hence the integrals on V and S n−1 make sense. In the next Lemma we study the biologically meaningful case of a network completely aligned in a single direction.
Proof. Since p ∈ ker L and since it is spatially homogeneous, equation (17) is satisfied. To study (18) we first compute the following integrals on S n−1 for all x ∈ R n .

Pointwise steady states.
In the preceding Lemmas we identified two simple homogeneous steady states. A full analysis of other steady states at this level is difficult, since the very weak formulation of measure valued solutions allows too many degrees of freedom. We rather specialize to the study of pointwise steady states as defined below. With pointwise steady states, we can combine the previous two Lemmas and design networks of aligned tissue with patches of uniform tissue. A schematic of the steady states which we construct here is shown in Figure 3.
4. For each test function Φ ∈ C(V ) and each Remark 1.
1. An immediate consequence of item 1. and 4. in this definition is that pointwise steady states satisfy for each test function φ ∈ W 1,1 (R n , C(V )). 2. Another immediate observation is that the homogeneous steady state from Lemma 4.1 and the completely aligned steady state from Lemma 4.3 are pointwise steady states.
In the following we classify pointwise steady states in R 2 . It turns out that the above definition allows for patchy steady states and for steady states of network type. Patchy steady states include patches of uniform tissue surrounded by areas of strictly aligned tissue, i.e. a combination of the above two types. Network type steady states arise if the areas of aligned tissue form a connected network of curves with intersections and branches. For network type steady states, we will classify possible intersections of network fibres. We are able to explicitly treat intersections of up to four directions and we find a general algebraic condition for networks with intersections of higher order.

4.2.
Patchy steady states. Assume a set of smooth closed curves σ i , i = 1, . . . , N separate R 2 into disjoint open sets Ω i , i = 1, . . . , k. Assume these curves σ i have finite length, no intersections but they might be closed. Assume p and q are uniform inside each patch with p i = 0 if |Ω i | = ∞. Since we are interested in pointwise steady states, we need to define (p(x), q(x)) for x ∈ i σ i . For each i = 1, . . . , N we denote the unit tangent vector at x ∈ σ i by γ i (x), where we will suppress the argument x whenever possible. We define and i ≥ 0.
Lemma 4.5. The weak steady state defined by (22) and (23) is a pointwise steady state.
Proof. (p(x), q(x)) are defined for all x ∈ R 2 and as shown in the proofs of Lemma 4.1 and 4.3 the conditions (19) and (20) are satisfied for all x ∈ R 2 . We only need to show that (p, q) as defined above is a weak steady state, i.e., we need to confirm condition (21). We find The first integral vanishes since ∇ x p(x, dv) = 0 in Ω i . The boundary integrals are zero, since we assumed that on σ i the fibre orientation is tangential, i.e. n · v = 0 for all v ∈ supp p(x, t), where n denotes the outer normal of σ i at x ∈ σ i . The above lemma allows for patches of uniform tissue surrounded by aligned tissue. These could also be called encapsulations, as seen for many tumours in tissue. A schematic of patchy steady states is given in Figure 3C. The steady states in Lemma 4.5 do, however, not allow for intersections of the curves σ i so that they become of network type. To obtain network steady states, we need to study possible intersections in more detail.
Lemma 4.6. Assume (p, q) is a weak steady state of (7) and at x it is of the form (24). It can only be a pointwise steady state, if the directions γ 1 and γ 2 have equal weight, i.e. if α = 1 2 . Proof. We only need to check condition (19) of Definition 4.4 at the intersection point x. For this choice of p and q we find and where we used the fact that |γ i γ i | = 1, i = 1, 2. To check condition (19), we need to test with a test function Ψ ∈ C(S 1 ): Hence to satisfy (19) for any test function, we need to satisfy Comparing the first and last term, we obtain the equation which is satisfied only if α = 1 2 . Notice that we assume |γ 1 γ 2 | = 1. For α = 1 2 we find B(p, q) = 1 2 and hence the condition (27) is satisfied. Next we study the general case where at a given point x ∈ R 2 we have an intersection of N -different directions γ 1 , . . . , γ N ∈ S 1 . We study N directions with equal weight: To decide if this intersection can be a pointwise steady state, we define a matrix of pairwise projections: Γ := (|γ i γ j |) i,j=1,...,N .
Theorem 4.7. Assume (p, q) is a weak steady state of (7) and at x it is of the form (28). It can only be a pointwise steady state, if the corresponding projection matrix Γ has an eigenvector (1, . . . , 1) T .
Proof. Again, we only need to check condition (19) of definition (4.4) at the intersection point. For the above choice of p and q we find and Applied to a test function Ψ ∈ C(S 1 ) we obtain and To satisfy condition (19) the right hand sides of (32) and (33) have to coincide for any test function. In particular we need to satisfy This condition implies |γ k γ i | for each l, k = 1, . . . , N.
It can be directly verified that condition (35) implies (34) and it also implies (32)=(33). Hence (35) is the limiting condition. This condition implies that the row-sums of the matrix Γ are all identical, and since Γ is a symmetric matrix, the column sums also have the same value. In other words (35) is equivalent with the statement that Γ has an eigenvector (1, . . . , 1) T . A schematic of a steady state with intersections is shown in Figure 3D.
Notice that a related matrix to Γ is well known in linear algebra: the Gram matrix G = (γ i γ j ) i,j plays a role in coordinate transformations and the square root of the Gram determinant is a measure for the volume element spanned by the vectors γ 1 , . . . , γ N . Example 4.8. As an example, we apply this general result to the two-directional case studied in Lemma 4.6. For two directions we have For three directions we obtain an interesting result: Assume (p, q) is a weak steady state of (7) and at x it is of the form (28) with N = 3. It can only be a pointwise steady state, if the three directions have equal angle, i.e. |γ 1 γ 2 | = |γ 2 γ 3 | = |γ 3 γ 1 |.
The classification of intersections of four directions is a bit more complex. (p, q) is a weak steady state of (7) and at x it is of the form (28) with N = 4. It can only be a pointwise steady state, if the pairwise equal angle condition (36) is satisfied.
Hence we obtain six unknowns and three equations, which will not give us such a complete solution as for three directions. We can, however, reduce the above condition to a set of pairwise equal angle conditions If all angles are π/2 then this condition is satisfied.

Unsymmetrical intersections.
In the numerical simulations shown in Figure 2 we observe intersections that are asymmetric in the sense that for a direction γ 1 the opposite direction −γ 1 is not seen. Indeed, at an intersection point x various fibres come together. This means that at this point x the network has a number of directions γ 1 , . . . , γ N , but the opposite directions are missing (it could of course happen by chance that γ j = −γ i for some i, j). Even though we assume that the distribution function q is symmetric almost everywhere, we find exceptional points at those intersection points. In this section we show that unsymmetrical intersections can arise as steady states in the framework developed here. Assume at a given point x ∈ R 2 we have an unsymmetrical intersection of Ndifferent directions γ 1 , . . . , γ N ∈ S 1 with equal weight: q(x) := 1 N (δ γ1 + · · · + δ γ N ) and p(x) =q(x).
To decide if unsymmetrical intersections can arise as pointwise steady states, we carry out the same computations as in the previous section. It turns out that the computations change only marginally and we omit the details here. For example formulas (32) and (33) use Ψ(γ i ) instead of Ψ(|γ i |). This implies that the conditions for their existence remain the same. We summarize: Theorem 4.9.
1. Assume (p, q) is a weak steady state of (7) and at x it is of the form (37). It can only be a pointwise steady state, if the corresponding projection matrix Γ has an eigenvector (1, . . . , 1) T . 2. If N = 3 then (p, q) can only be a pointwise steady state, if the three directions have equal angle, i.e. |γ 1 γ 2 | = |γ 2 γ 3 | = |γ 3 γ 1 |. 3. If N = 4 it can only be a pointwise steady state, if the pairwise equal angle condition (36) is satisfied.
1. Indeed, unsymmetrical intersections of three directions with angles of 120 • seem to be typical building blocks for the network shown in Figure  2. 2. Although other intersections (symmetric or asymmetric) do exist theoretically, they are rarely seen in simulations. This raises the question of stability of these steady states. We defer this question to future studies.

4.5.
Other steady states. We consider more general steady states where the cell distribution is a multiple of the lifted fibre distribution, that is where ∈ L ∞ (R n ) is the density of cells (or even ∈ L 1 ∩ L ∞ (R n )). The minimal condition for such a pair to be a steady state is in particular, the left hand side is actually independent of θ. Because of the linearity of the operators Λ and B, this condition becomes (x)Λ(q(x))(θ) = (x)B(q(x), q(x)).

Wherever = 0 this condition can be stated as
for almost all x ∈ R n . Now let us try the following ansatz where 0 ≤ f (x) ≤ 1 and Σ is the normalized Haar measure on S n−1 . Notice that even the predominant direction γ may depend on x at this point. A calculation gives The last term in the second line is independent of θ because of the rotational invariance of Σ. Hence we have Observe that (x)q(x) = (x) and the right hand side of the first equation of (7) vanishes. Hence the condition for is (x)∇q(x) +q(x)∇ (x) = 0, for almost all x ∈ R n , which can be written as In summary, if we have found q that satisfies (38), then can be determined from the differential equation (39).

5.
Analysis for given fibre distribution. In this section we assume that cells do not remodel the fibre network and that q(x, t) is a given distribution. The pequation of system (7) has a simple structure for given q. For classical solutions we can use the method of characteristics to find an explicit solution. For given v ∈ V , the characteristic equation is d dt x(t) = v. Hence the characteristic through x 0 ∈ R n is given by x(t) = x 0 + vt. We can write the first equation of (7) as follows We evaluate equation (40) at V and obtain d dtp Applying this to equation (41), we have This is an equality in the Banach space B(V ) and the term p 0 (x − t dv) on the right hand side has to be interpreted as the v-shifted measure defined in equation (9). The same notation applies toq. We evaluate this measure p(x, t) at V , i.e., we computep(x, t) = p(x, t, V ), and obtain which is an equality between real numbers. The measureq is non-negative and for fixed w ∈ V we haveq(x − w(t − s), s)(V ) = 1, and Thus we get (1 − K(x, t))p(x, t) = e −µt p 0 (x − V t, V ).
If K(x, t) = 1, then we can solve forp as Thenp can be used in (42) to find an explicit solution Notice that this solution only depends on the initial condition p 0 and on the fibre distribution q. To clarify, equation (45) is again an equality in B(V ) and the right hand side is of the type "measure + number · measure". Equation (43) simplifies drastically in the special case of constant fibre distribution q(x, t) = q. In this case, it follows that Equation (46) deserves some interpretation.p is the mass density of particles of all velocities at point (x, t), whereas p 0 (x−V t, V ) integrates the initial condition over the domain of dependence of the point (x, t), the set {x − tv : v ∈ V }. The velocity distribution at (x, t) arises by following all characteristics through (x, t) backwards (see Figure 4). We call (46) a generalized Huygens principle. The solution p(x, t) from (45) can then be written entirely in terms of the initial condition Using equation (46), the explicit solution can also be written as Hence the solution is a convex combination of the initial condition p 0 and the current amount of cellsp redistributed with respect to the "controlling" distributionq. We will use this observation in the next section to rigorously prove convergence of the parabolic limit.
x,t 1 x t x,t 2 It is interesting to understand the biological meaning of these explicit solutions. Equation (48) tells us that eventually cells will completely align to the given network structure. This has similarities with glioma cell invasion along white matter tracks of the brain [17]. Equation (45) has a similar interpretation. In contrast to (48), here the fibre distribution varies in time and space. The integral terms in (45) and (43) denote a temporal average over the history of the fibre orientation, where the influence of the history is exponentially damped. Then (45) can be understood as cells that try to align with the tissue while the tissue is changing. We see a biological analogy with wound healing, where fibroblasts constantly modify the collagen network while immune cells move through this fibre scaffold and heal the wound [7]. Notice that it is not our intention here to model brain tumours or wound healing. These examples are only used as analogies and thought experiments. It might be useful to make our model available to these processes in the future.
5.1. The parabolic limit problem. As shown in [14], we can formally derive a diffusion limit equation from equation (7) under suitable scaling of space and time. Letx andv denote reference length, and speed, respectively, with the dimensionless quantity ε =v µx being small. We introduce rescaled variables as follows
(P ε ) Simultaneously, we consider the limit problem with || ( · , 0)|| ∞ < ∞ and where the diffusion tensor is given by The formal derivation of the limit problem and the diffusion tensor in equation (49) was carried out in [14, section 4], see in particular equations (29) and (41) in that paper. We therefore omit these calculations here. Notice that D[q] can be written as the scaled variance-covariance matrix V(q) of q, and we have used equation (3).
By the symmetry of q, we have We divide equation (53) by ε 2 and let ε → 0 and obtain ∂ ∂t where ∂p ε ∂t ∂ ∂t in the distributional sense. Using the representation (52) we obtain Hence satisfies the limit equation (P 0 ).
6. Discussion. In this paper we consider mathematical properties of a model (1), or equivalently (7), that describes mesenchymal cell movement in tissues. The model was developed in [14] and has been analyzed from various angles in recent papers, [21,4,5,25]. Through the previous analysis it became evident that a solution framework is needed which allows for measure valued solutions. Here we develop such a framework and prove global existence of solutions in X. We have used semigroup methods, since they provide a dynamical systems point of view, and we can use this framework for linear stability analysis in future work. Alternative methods to show existence include energy methods as developed by DiPerna and Lions [9]. We were able to find non-trivial measure-valued steady states, which correspond to homogeneous distributions, or to aligned tissue, or to patches of uniform tissue with a network separating these patches. We found that pointwise steady states show network properties as observed numerically. We also found that, although the system has been formulated symmetrically, we can have unsymmetrical intersection points. This confirms the interpretation in [14], where it was suggested that a network made from undirected fibres can have characteristics of a directed network. The complete identification of steady states of (7) is an interesting open question. Furthermore it would be interesting to see whether solutions of (7) converge to steady states or to traveling wave solutions as t → ∞. The existence result in X opens the door to a rigorous linear stability analysis of steady states. This endeavor is left to future work. The convergence to the parabolic limit is a standard feature of kinetic models and it has been studied in many publications (see references in the text). Our approach here extends known results to measure-valued solutions. Furthermore, we formulate an explicit solution which shows that the solution basically is a convex combination of initial data and its velocity-mean-value. The mean value then is close to the parabolic limit. We also give an argument that the rigorous convergence to a diffusion limit might only work for constant tissue.
Here we did not discuss the biological modelling of (7). We would like, however, to discuss the biological assumptions and propose various extensions, which could lead to more realistic models.
One possibility is to introduce birth and death processes for the cells into the model. For example, it is known that growth factors can be bound to the fibres which promote proliferation. Also, harmful substances can be found in the fibre network, possibly killing the cell. To model these effects, a term G(p) would be added to the right-hand side of the first equation of (7). If G(p) satisfies certain growth bounds, the global existence of solutions will continue to hold.
A second possibility is to allow diffusion of p with respect to both x and v. Cells are likely to undergo some random walk and may also change their velocity randomly (a perfect alignment of cell velocities will disperse). Diffusion with respect to the x variable is easily modeled by adding a term of the type −D x ∆ x p to the left-hand side of the first equation. Also, diffusion in the velocity can be modelled through an additional diffusion term of the form −D v ∆ v p (see also Dickinison [8] for chemotaxis). For these cases we expect a smoothing property of the linear semigroup and the totally aligned steady states will no longer exist.
Another possibility to expand and make the model more realistic would be to give the fibres some elasticity and to let the fibres be moved by the cells. Also, the cells should chose their new speed not randomly, but according to some "stiffness" of the neighborhood they are currently in. For example, a cell that has to cut a lot of fibres in its way should slow down, while a cell that is aligned well with the network can gain speed. Obviously, these are intuitive ideas, and would have to be supported by biological evidence.
In model (7) we implicitly assume that the protease is released locally at the leading edge of the cell. In the literature, however, various protease cutting mechanisms are discussed [11] and more detail of the cutting can be included into the model (see also Painter [21]). This might necessitate to explicitly model the protease as a third variable through its own reaction-diffusion equation.
A consideration of the length scales of the fibres relative to the size of the moving cells might also give valuable input into the appropriate modelling assumptions.
Finally, we have studied an unbounded domain R n to avoid boundary conditions. To formulate the correct boundary conditions for model (7) is not trivial. A commonly observed effect seen in tissue is that a tumour is encapsulated by a dense fibre network. In that case the fibres at the boundary will be aligned tangentially to the boundary and should trap moving cells inside the domain. The encapsulations can be understood as patchy steady states, as described here. A careful analysis of other boundary conditions and its implications on existence and steady states is left to future work. For the simulations in Figure 2, K. Painter used periodic boundary conditions on a square domain, i.e. a flat torus.