Anisotropic Network Patterns in Kinetic and Diffusive Chemotaxis Models

: For this paper, we are interested in network formation of endothelial cells. Randomly distributed endothelial cells converge together to create a vascular system. To develop a mathematical model, we make assumptions on individual cell movement, leading to a velocity jump model with chemotaxis. We use scaling arguments to derive an anisotropic chemotaxis model on the population level. For this macroscopic model, we develop a new numerical solver and investigate network-type pattern formation. Our model is able to reproduce experiments on network formation by Serini et al. Moreover, to our surprise, we found new spatial criss-cross patterns due to competing cues, one direction given by tissue anisotropy versus a different direction due to chemotaxis. A full analysis of these new patterns is left for future work.


Introduction
Network formation is an important process in many biological tissues. For example, in angiogenesis, blood vessels sprout and grow to form functioning vascular networks [1], or lung tissue, which forms a complex network of bronchioli. Another example is the developing brain, where glial cells grow cellular protrusions, which eventually mature into a synaptic neuronal network [2]. In addition, on the pathological side, we see network formation, for example, in brain cancer growth [3,4]. Thus, a good understanding of network formation can help to identify healthy and pathological processes and provide information for control through treatments. Mathematical modeling has played an important role in physiological processes for many years, and the modeling of network formation is no different. The mathematical approaches are manifold, involving agent-based models, stochastic processes, and continuum models [3,5]. In most of these models, a leader cell (or tip cell) invades new areas, leaving a trail behind, which follower cells can identify and follow. Eventually, the trail forms connections and matures into a network.
A different mechanism for network formation acts in vasculogenesis [6][7][8]. The spatial interaction of moving endothelial cells leads to a spontaneous network formation, in a process that is not based on a trail following mechanism but rather arises as macroscopic pattern formation in a chemotaxis system. Serini et al. [6,7] performed experiments and developed a mathematical chemotaxis model to describe this process. Serini's et al. model is based on a macroscopic population-level description, which we discuss in a bit more detail in Section 2. Here, in this paper, we generalize Serini's et al. approach and consider a chemotaxis-transport equation model, which is based on detailed microscopic cell movement patterns often referred to as "run and tumble" [9][10][11][12][13].
There are several ways to include chemotaxis, the active orientation of individuals along chemical gradients, into continuum models [14]. In the case of transport equations of diffusion type, chemotaxis was included by Alt [10] and analyzed in detail by Hillen and Othmer [12,13]. Chalub et al. [15], and Bellomo and Winkler [16] used the kinetic theory for active particles [17] to derive chemotaxis models. Chalub et al. [15] develop a complete existence and solution theory for these models, and their methods can be applied to our model below with minimal modifications. The details of the existence results have been worked out in great detail in Thiessen [18]. Here, we follow the idea of Kumar and Surulescu [19], where we derive a fully anisotropic chemotaxis model from an anisotropic transport equation framework. The chemotaxis-transport approach has been recently applied by Kumar et al. [19,20] to describe tumor-palisades. Tumor palisades arise as microscopic patterns in certain brain cancers, where cancer cells arrange themselves in a ring-like structure. In Kumar et al. [19], it is shown through mathematical modeling that these structures are likely a result of increased acidity in the tumor center.
In the general chemotaxis-transport framework, we repeat a similar process to Kumar et al. and once again retrieve anisotropic diffusion. Using parabolic scaling methods, we find an anisotropic chemotaxis model, which shows network forming properties, but not only that. The anisotropic chemotaxis model balances two causes of anisotropy, the directionality of the underlying tissue and the orientation of the chemotactic gradient. These two effects compete, leading to interesting criss-cross patterns, which have never been observed before in a mathematical model. Accurate numerical simulation of the anisotropic chemotaxis model is a challenge, and we develop an efficient finite volume finite difference solver (FVFD). The numerical solver was inspired by the work of Chertock et al. [21], where a second-order finite volume finite difference solver for a closely related Keller-Segel system was developed.
Our results give new insights into the interplay of two-directional cues on moving populations. The directional effects do not arise as a naive linear combination of each individual effect; rather, the anisotropy and chemotaxis sensing can interact in complicated ways. Further analysis of mutliple directional cues in a transport equation framework is currently been investigated by Conte and Loy [22].
The paper is organized as follows. In Section 2, we revisit the model of Serini et al. [6] and we give some background material on transport equations in biology. In Section 3, we develop the transport-chemotaxis model, and we perform the parabolic scaling to derive the anisotropic chemotaxis model (34). In Section 5, we develop an efficient numerical solver, and we apply it to some relevant examples. Here, we show the occurrence of the new criss-cross patterns. We close with a Discussion Section 6. In 2003, Serini et al. [6] performed experiments on vascular endothelial cells and developed a continuum model to explain vascularization. Endothelial cells were cultured on a matrigel surface, and over time network forming patterns were observed. These patterns are controlled by the total cell density and by the strength of their chemotactic response. Small cell densities lead to local cell alignment, intermediate cell densities lead to cell network formation, while large cell densities lead to broadly aligned "highway" networks.
From the experiments in Serini et al. [6], we can identify time and space scales, which will inform the parabolic scaling, which we do later. From the videos and images provided, we find a microscopic time scale of 1 min, where cells change direction about once every 4 min. During an observation of 1 min, they travel a distance of maximal 0.05 mm. The macroscopic scale of observation is about 5-9 h on a domain of about 1 mm 2 . If we introduce a small parameter ε = 0.05, then we find the ratios microscopic space macroscopic space The mathematical model of Serini et al. [6] is a system of three partial differential equations for the balance of cell density, macroscopic cell population momentum and the chemoattractant. The cell population is described by a continuous density ρ(x, t) and population movement is described by a macroscopic velocity v(x, t). Gradients of a chemoattractant S(x, t) act as accelerating force on the population velocity, while the chemoattractant S is released by cells, diffuses through the environment and degrades over time. The resulting equations are where β is a constant denoting the strength of chemoattractant sensing, D is the diffusion constant of the signal, α is the rate of chemoattractant production, and ϑ is the decay rate of the chemoattractant. Finally, the term p denotes a pressure term, which is used to regularize the solutions and avoid unrealistic density blow-up. The model reproduces the observed network patterns very well, confirming that network patterns only formed with a cell density of at least 100 cells/mm 2 .
Here, we extend the Serini et al. [6] model (2) in two ways. Firstly, instead of using a macroscopic population-level description for mass and momentum balance, we go back to the behavior of individual cells. We develop a transport equation model that includes movement characteristics, such as cell velocity, turning frequency, and choice of preferred directions. Secondly, our formulation allows for a very natural inclusion of tissue anisotropies that arise from aligned tissue structures, such as collagen fibers, nerve fiber bundles, etc. Coupled with chemotaxis, we obtain an exciting interaction between cell movement, tissue anisotropy, and chemotaxis. The main focus of our paper is to better understand this interaction.

Transport Equations
Transport equations have become a powerful tool for spatial modeling in biology. They are particularly useful in cases where individual particles, such as cells or animals, can be followed and characteristics about their path, such as speed, direction, and turning rate, can be measured. In the context of this paper, we are concerned about two forms of cell alignment. On the one hand, cells align with linear features of the underlying tissue, such as blood vessels or nerve fibers or collagen fibers. In addition, cells react to a chemoattractant and orient themselves in the direction of increased concentration. A transport equation framework is an ideal environment to develop models for oriented movement.
Transport equations describe the time evolution of a particle density P(t, x, v) at time t ≥ 0, location x ∈ Ω ⊂ R n and velocity v ∈ V ⊂ R n . We consider a periodic bounded domain Ω, or the whole space Ω = R n . We introduce a set V to denote the possible velocities of the particles and we usually take V = [s 1 , The time evolution of P(t, x, v) is described by the transport equation where the index t denotes the partial time derivative, and L is the turning operator, which describes the specific directional changes of the particles [9,10,12,23]. L describes how particles chose a new direction, how strongly they are influenced by the underlying directional information, and how persistent their movement is. Typically, L is defined via an integral operator representation: where the first term on the right-hand side gives the rate at which particles switch away from velocity v, and the second term denotes the switching into velocity v from all other velocities v . The spatially dependent parameter µ(x) is the turning rate, while 1/µ(x) is the mean run time at location x. The traditional literature has focused on the case of constant turning rate [12,23], while direction-dependent turning rates have recently been studied in Reference [24]. The kernel T(x, v, v ) denotes the probability density of switching velocity from v to v, given that a turn occurs at location x. The properties of the turning kernel T are the key to the following analysis, and many different choices are possible. Roughly, we can classify the most important cases into three categories, where we will focus on the last one in this manuscript.

1.
Boltzmann equations: Transport equations, such as (3), have a long history in thermodynamics and physics [25]. In that case, the turning kernel L is a non-linear collision kernel between gas particles (or other particles) and the set of velocities V is unbounded. On the other hand, a linear kernel, such as (4), arises as linearization on a so-called Maxwellian distribution [25]. We see later that an important difference to biological applications is the fact that the Boltzmann equation conserves mass, momentum, and energy, i.e., the kernel (null space) of the interaction operator L is five-dimensional. In our applications, we have only one conserved quantity: mass, and only in cases where there is no particle birth or death.

2.
Diffusive transport equations: The diffusive transport equation is characterized by a turning operator L, whose kernel consists of functions that are constant in velocity. This means, in equilibrium, there are no preferred directions, and on a long time scale, the dynamic becomes diffusion-like. In References [12,13], we developed a complete theory of diffusive transport equations, and we have shown how chemotaxis models can arise within this theory.

3.
Anisotropic transport equations: In this case, particles have strong orientational guidance. This guidance arises either through a preferred equilibrium distribution, called the Maxwellian, or it arises from external directional cues. The models in this class fall within the kinetic theory for active particles (KTAP) [16,17]. Mathematically, this case is characterized via a non-trivial one-dimensional kernel of the turning operator L. The case of an external directional guidance is of interest here, and a whole theory for external directional cues was developed in References [26][27][28]. We present this case in more detail.
In the anisotropic case, it is assumed that the turning operator does not depend on the incoming velocity v , i.e., where q satisfies q ≥ 0. This assumption seems restrictive since many species would show some form of persistence in the movement direction, and it is entirely possible to include such a dependence into an anisotropic movement, as well. However, it has been shown in many applications that the simplifying assumption (5) is extremely useful in the modeling process, and it can be justified in many cases (certainly for cancer invasions [29] and for sea turtle navigation [30]). The anisotropic framework was developed in Reference [26] and extended in Reference [27]. We consider a given environment, where at each location in space a distribution of preferred directionsq(x, θ) is given, where θ ∈ S n−1 is a unit vector. We assumẽ We assume that particles that are changing direction do so by using the underlying network structure, i.e., we assume that q(x, v) ∼q(x,v), wherev = v/||v|| denotes the corresponding unit vector. Sinceq is a probability distribution on S n−1 , and q a probability distribution on the set of velocities V, we need to normalize appropriately: This definition motivates us to abuse notation and call q the directional distribution of the underlying environment. In doing so, we need to remember that q is a rescaled version ofq, which is really the directional distribution of the underlying environment.
For this choice of turning kernel (5), Equation (4) simplifies to To summarize, the anisotropic transport equation has the form For the following analysis, it is useful to consider two statistical quantities of q, the expectation and the variance-covariance matrix of q(x, v) on V: Let us discuss a few examples of directional distributions q.
For these examples, we assume that the particle speed is constant s, i.e., V = sS n−1 . The extension to V = [s 1 , s 2 ] × S n−1 is straightforward, but it introduces some extra parameters that need to be carried through the integrals.

1.
Uniform distribution: The simplest case arises if there is no directional cue at all. In this case,q This uniform case is also included in the class of diffusive transport equations mentioned above. In this case, where I denotes the identity matrix.

2.
Strictly-aligned tissue: If all fibers are pointing in the same direction γ(x) ∈ S n−1 , In this case,

von-Mises distribution:
The von-Mises distribution is the analog of the normal distribution on the unit sphere S n−1 . Strictly speaking, it is called von-Mises distribution only on S 1 . In higher dimensions, it is called the Fisher distribution. However, the functional form is the same as in 2-D [31]. In 2-D, the von-Mises distribution is where k(x) ≥ 0 is the parameter of concentration, γ(x) ∈ S 1 are given preferred directions, and I 0 (k(x)) is the modified Bessel function of first kind of order 0. As k → 0, the von-Mises distribution becomes uniform (case 1), and, for k → ∞, the von-Mises distribution becomes singular (case 2). The expectation and variance of the two-dimensional von-Mises distribution are [31] E q (x) = s where I 0 , I 1 , I 2 denote the modified Bessel functions of the first kind of order 1, 2, and 3, respectively. 4.
Bimodal distributions: In many cases, q is symmetric, and guidance up a fiber or down a fiber has the same probability, with expectation and variance [31] E q (x) = 0 , Note that, instead of using a dyadic product γγ T , some authors prefer to use a tensor product γ ⊗ γ. We have no preference, but we chose to use only one of these notations.

Chemotaxis-Transport Equations
As discussed in the Introduction, chemotaxis effects have been included into the kinetic framework by several authors [10,12,13,[15][16][17][18][19]. We now extend the model (6) to include the effect of chemotaxis. For simplicity, we are considering the spatial domain to be the torus, which we will denote as Ω = T n . The turning operator only involves cells changing velocity and, therefore, should not affect the number of cells. Assuming that the number of cells remains the same, we gain the constraint on the turning operators The corresponding conservation law can be obtained by integrating (3) over the velocity domain where we use the density of cells and the first moment as
To include chemotaxis into the transport equation framework, we assume that the turning kernel T depends on the cell density and the gradient of the chemotactic signal S. T := T[ρ, S](x, v). The movement of the chemical signal S(x, t) is traditionally defined as a diffusion equation. However, depending on the time scale at which the chemoattractant diffuses, one can consider that chemoattractant is at equilibrium, then the resulting concentration is given by a Poisson type equation. Under this fast diffusion assumption, powerful results on existence and uniqueness have been obtained [32]. We take the traditional route and describe the motion as a diffusion equation, where the chemoattractant production is proportional to the cell density, and the chemoattractant degrades over time. These assumptions give us the diffusion equation Here, S := S(x, t) is the concentration of the chemoattractant, D s is the diffusion constant associated with how the chemoattractant moves through the environment, α is the rate of production of the chemoattractant by the cells, and ϑ is the rate of degradation of the chemoattractant.
For our purposes, we specify the turning operator further. Inspired by Reference [13], we split the two guiding effects into the environmental contribution q and chemotactic effects as Here, we allow the environmental cues to depend on the cell density q[ρ](x, v). The rate b[S](x, v) describes the chemotactic strength, based on the concentration of the signal S. In most examples, b is a constant. The factor (1 − ρ) is a volume filling term, and it allows chemotactic movement only into areas where there is space available. The variable v denotes the outgoing velocity, and ∇S is the gradient of S. To ensure that T is a probability distribution, we need V vb[S](x, v)dv = 0, which we generalize into the following assumptions: The examples we will be looking at is the von-Mises distributions (8) or (11) for q and where β is the general sensing strength, and A(x) is a matrix that describes the directional dependence on sensing the chemoattractant. Note that the above form of b(v) is a general second order expansion of a function of v with vanishing odd moments. We can take advantage of the structure we have provided, and the turning operator Inputting this turning operator into (3), we obtain the main chemotaxis-transport system of this paper: We now have a mesoscopic set of equations; in the coming chapters, we will take the parabolic scaling of these equations to get a macroscopic system and explore numerics to get insight into pattern formation.

Parabolic Scaling
We have seen above that the transport Equation (20) acts on microscopic scales. For example, in the case of Serini's experiments [6], we have a microscopic time scale of minutes and a microscopic length scale of about 0.05 mm, while the macroscopic scales are several hours and millimeter. This situation allows us to use parabolic scaling to derive a model that lives on the macroscopic scale. The parabolic scaling of transport equations is a welldocumented method, and explorations can be found in References [10,12,15,16,19,23], for example. Similar to the scaling (1), we assume there is a small parameter ε > 0 such that space and time are scaled as Then, the above transport Equation (20) becomes To analyze the scaled Equation (22), we take the scaled coordinates (ξ, τ) and make a regular expansion in ε, called a Hilbert expansion [33]: A further assumption is that all of the mass is contained in the first order Substituting the Hilbert expansions into the rescaled transport Equation (22) and comparing terms of equal order yields a countable number of equations at different orders ε that all must vanish independently. The zeroth-order equation is given by The next order of ε gives Solving this equation in terms of P 1 leads to To get a closed system for the first term in the Hilbert expansion P 0 , we look at the ε 2 terms Integrating over V, the right-hand side vanishes, and substituting P 1 from (27) and P 0 from (25) yields an equation in terms of ρ To simplify this equation, we temporarily use index notation, where, according to Einsteins summation convention, we sum over repeated indices.
Using the previously defined expectation and variance (7), we can write the above equation in a more compact form with the identity We use the same definition of expectation and variance for other functions; for example, Combined with the above chemotaxis Equation (21), we obtain the macroscopic anisotropic chemotaxis system whereα = ε −2 α, andθ = ε −2 ϑ. In other words, on the microscopic scale, the rates α and ϑ are small, and, on the macroscopic scale, they are of order one. Since we are solving the chemotaxis model only on the macroscopic scale, we remove the overline-tilde's and simply keep α and ϑ.
The system (34) is a generalization of the typical Keller-Segel equations [34] in the sense that we include anisotropic diffusion V q + E q E T q , and anisotropic chemotaxis V b . The effect of the V q and V b terms are nonuniform diffusion and mixing of the ∇S into other directions (chemotactic mixing). Of course, the above scenarios are dependent on the given distributions q[ρ](x, v) and b[S](x, v). For instance, if q and b are uniform distributions in v, then the above model is the fully isotropic Keller-Segel system [14].

Examples
To explore the effects of anisotropy, we consider the examples that we mentioned above. We use a bimodal von-Mises distribution for q as in (11) and a second-order expansion for b from (18). The moments of the bimodal von Mises distribution are given in (12). For example, when the fiber direction is in the top-right diagonal direc- The first term is isotropic diffusion while the second term causes diffusion to be increased along the directions ( 1 I 0 (k) → 0 as k → 0, making the diffusion tensor isotropic, and the case I 2 (k) I 0 (k) → 1 as k → ∞, making the diffusion tensor to be maximally anisotropic.
To compute the second moment for b from (18), we use index notation and summation convention, and we ignore the explicit x dependence: We move the constants out of the integrals and define a general mean velocity tensor to bev Then, we can rewrite and, using the Lemma 2.2 from Hillen [35], we explicitly compute thev's from (36) as where the set of pairs of indices out of (i 1 , . . . , i 4 ) is defined as We compute the sum explicitly since there is only 12 combinations, and, due to the symmetry of the Kronecker delta, there are only 3 unique terms: Using the properties of the Kronecker delta, we can convert back to matrix notation As a result, we have anisotropic components caused by the off-diagonal pieces of A. If A has nonzero off-diagonal components, then V b has a mixing effect on the chemotactic velocity u(x, t) := V b ∇S.
To illustrate the chemotactic mixing in two dimensions, we consider the velocity in the x direction is dependent on the chemotactic gradient in the y direction, and vice versa. In summary, the new terms create non-uniform diffusion and mixing of influence of the chemotactic gradient. The effect of the chemotactic mixing will be shown through numerics in the next section.

Numerics
In this section, we numerically explore the effect of anisotropy on chemotaxis. For this exploration, we return to the example of the endothelial cells randomly distributed to form a vascular network. We describe this motion using the parabolic limit (34). For the domain, we consider the two-dimensional square [0, L] × [0, L] with a periodic boundary condition. We discretize this square into smaller squares with side lengths of h and index these squares by the indices j, k, representing the x-coordinate and y-coordinate, respectively. Both indices run from [1, N], where N = L h .

Hybrid FVFD Scheme
For the numerical scheme, we will closely follow the method proposed by Chertock et al. [21]. Chertock et al. derived a second-order hybrid finite volume finite difference (FVFD) scheme for the classical Keller Segel model They were able to show positivity preservation for this scheme based on reasonable CFL conditions, making it ideal for our problem. The numerical scheme is based on the Godunov scheme for conservation laws [36] to ensure the conservation of particle mass. Since our system has added complexities through the anisotropic diffusion terms, we need to take special care identifying the appropriate upwind terms in the numerical flux. For this reason, we will go through the scheme in detail to exhibit these differences. To begin, we have to write (34) in the form of fluxes along the axes; for this purpose, we write (34) in two dimensions: On the above mesh, the goal is to write the system as a semi-discrete hybrid FVFD scheme, which takes the form of where the averaged cell density isρ j,k ≈ 1 ρ(x, y, t)dxdy, and S j,k is the point value of the chemoattractant. F j,k are the numerical fluxes in the x and y directions, and ∆ j,k is the discrete Laplacian.
We discretize the ρ equations with a finite volume method since finite volume methods are conservative. The conservative property of this scheme can be seen by looking at the conservation law where u is the state of the solution, and F(u) is the corresponding flux vector. Now, we divide the domain of this problem into finite volumes V i , and integrating over one such cell yields On integrating the first term to get the volume average (like the approximation ρ j,k ≈ 1 h 2 I j,k ρ(x, y, t)dxdy) and applying the divergence theorem to the second yields where S i is the surface of the volume of V i , and n is the outward normal. Now, since the change of cell averages is determined by the edge fluxes, any gain of one cell is a loss for another. Dividing by V i yieldsū In the case of two dimensions, we just need to keep track of the fluxes at the four edges surrounding a point. In our case, we subdivide our domain into squares with side length h centered at each grid point, and then the change in state cell average is where j + 1 2 , k, j − 1 2 , k, j, k + 1 2 , and j, k − 1 2 represent point values on the dual grid of spacing 1 2 h (see Figure 1). Based on Chertock et al.'s method [21], we write the numerical fluxes as follows: where u and v are the chemotactic velocities, and the point values ρ j+ 1 2 ,k and ρ j,k+ 1 2 are computed in a upwind manner as The one-sided point values, at the interfaces ρ E j,k , ρ W j+1,k , ρ N j,k , and ρ S j,k+1 , are computed using the second order piecewise linear reconstruction To ensure that the above point values are second-order and non-negative, slopes are computed adaptively: with the flux limiter minmod The positivity of the reconstructed cell density point values is guaranteed by the positivity preserving generalized minmod limiter [37][38][39][40], with the assumption that the underlying cell averages are positive. Now, the scaled cell density derivatives are computed via central differences Note that, in addition to Chertock's scheme [21], here, we have anisotropic flux terms which need special attention. We compute those as a finite difference over a double-spaced interval 2h since this allows us to use point values on grid points that are already computed (see Figure 1). Note that, in this formulation, we do not need values on diagonal points of the form (j + 1 2 , k + 1 2 ), making the scheme a bit easier. The chemotactic derivatives (S x ) j+ 1 2 ,k , (S y ) j+ 1 2 ,k ,(S x ) j,k+ 1 2 , and (S y ) j,k+ 1 2 are computed similarly.
For the discrete Laplacian, we use the standard five-point stencil [41] to obtain a second-order approximation Thus, we have derived a second-order semi-discrete method for (34): To evolve this semi discrete scheme through time, we chose to use a second order adaptive Runge-Kutta scheme [42].

Parameters
The system has a couple of parameters we have to determine namely, α the production rate of the chemoattractant, ϑ the degradation rate, µ the turning rate, D s the diffusion rate of the chemoattractant, and the forms of the distributions ρ(x, 0), q(x, v), and b(x, v). With the purpose of looking at vascular assembly, we can find the constant parameters, ϑ −1 = 3600 s, α = 1 from Reference [6], D s = 10 −7 cm 2 /s from Reference [43], and µ = 1785.71/s from Reference [44]. We take the initial condition to be a linear combination of Gaussian's where {(x l , y l )} W l=1 is a sequence of random numbers drawn from the uniform distribution on [0, N] 2 . This initial condition represents a random distribution of cells on a petri dish, which have a radius of σ. From available data of molecular radii [45,46], we can estimate the radius of a cell to be σ = 0.003 cm. With the general parameters defined, we can move onto the numerical experiments.

Scheme Test
In this numerical experiment, we attempt to show that our scheme's significant features are present for varying grid sizes. To show the invariance with respect to grid sizes, we take a single Gaussian in the center as the initial condition. As for the fiber and sensing directions, we use the von Mises distribution for q (11) and the quadratic form for b (18), as discussed earlier.
with diagonal anisotropy Calculating the variances (35) and (40), we find V q = 1 2 I + 1 2 where the max velocity s = 3 20π cm s , and the concentration parameter is k = 10. In Figure 2, we show the initial condition in Figure 2a and two simulations at time = 0.5 with different grid sizes. The grid size in Figure 2b is 100 × 100, and, in Figure 2c, it is 200 × 200. The two computations are qualitatively identical. These computations show the two main properties of our system. First, the cell density diffuses anisotropically along the main diagonal direction, which results from our choice of V q . Secondly, the peak splits into two aggregations near the center, which results from the chemotactic mixing V b at highly dense regions. In the following numerical experiments, we will further explore these aspects of the system and show that they lead to network formation.

Vascular Network Formation
In this numerical experiment, we are interested in replicating previous models' results in looking at vascular network formation as considered in Serini et al. [6,7]. In their experiments, they took the randomly distributed cells, with b being constant. For the distribution q, the choice is a little more complicated. In Serini et al. [6,47], they propose a pressure term p(ρ) = ρ 3 to deal with cell compression. They derived their models through the moment closure taking the form where Ep is the first moment of the pressure p. We match the pressure term with the strength of our anisotropic diffusion by setting V q = ρ 3 I. One of the critical findings of References [6,47] was that the formation of networks arises only for cell densities larger or equal to a critical value. For the above parameters, this value is 100 cells/mm 2 (a) To try to recreate the experiment in Reference [6], we evolve the Equation (34) from the initial condition (51) with the parameters ϑ −1 = 3600 s, α = 1, D s = 10 −7 cm/s 2 , µ = 1785.71/s; in addition, we take V q = ρ 3 I, V b = I.
In Figure 3, we show simulations of network formation for increasing total cell density. The top row shows the random initial condition, while the second row shows a snap-shot at time t = 1.5. In all three cases, localized structures are forming. For cell densities of 100 cells/ mm 2 , the patterns still appear as local aggregates; however, for larger cell densities of 200 and 400, larger-scale structures arise that resemble networks. Notably, we were able to obtain qualitatively similar results to those in References [6,47], despite using a different class of equations. One explanation lies in a moment closure approach that can be performed in the model of (53) and (54). Doing the moment closure in their model will lead to our macroscopic Equation (34) in the isotropic case (i.e., q and b are uniform in v). Details on moment closure methods are discussed elsewhere and are not part of our study here (see References [28,47]).

Anisotropic Diffusion
For this numerical experiment, we are interested in looking at the effects of anisotropic diffusion. We evolve the Equation (34) from the initial condition (51) with the parameters (55), and to explore anisotropic diffusion, we take ), which means that the variance matrices (35) and (40) take the form Cell density cells/mm 2 100 200 400 ) severely changes the system dynamics. There is still network formation, but now it is biased into the direction of anisotropy. As we increase the concentration parameter k, there is a tighter and tighter spread along the diagonal. We can also localize anisotropic effects to certain regions of the domain by taking the concentration parameter to be a function of x. For example, in the third column in Figure 4, we keep anisotropy in the center of the domain while removing it towards the boundary, i.e., We clearly see the effect of anisotropy in the domain center, while the regions near the boundaries behave isotropic.

Chemotactic Mixing
In this numerical experiment, we evolve the Equation (34) from the initial condition (51) with the parameters (55), and we are going to explore mixing of the chemotactic velocities through an anisotropy term in b. As seen at the end of Section 4, chemotactic mixing is where the chemotactic gradient in one direction affects the movement in all directions. For this effect to occur, we require non-zero off-diagonal components to V b ; therefore, we take Concentration parameter From Figure 5, we can see that the effect of the chemotaxis mixing is comparable to the effect of non-uniform diffusion as we have seen in the previous numerical experiment. The difference is that the chemotactic mixing seems to affect the diffusion along the diagonal (−1, 1), but the chemotaxis term usually acts as a deterrent, reducing movement away from high-density regions. Therefore, we hypothesize that the chemotaxis term is not promoting movement along the (−1, 1) direction, but discouraging diffusion along (1, 1), while leaving diffusion along (−1, 1) unimpeded, giving the structures in Figure 5.

Anisotropic Diffusion vs. Chemotactic Mixing
Now, we wish to examine the interaction between the fiber distributions q and the chemotactic velocity mixing. For this, we consider the q in (57) with the concentration parameter k = 5, and b in (18). According to the previous numerical experiments, these choices of distributions should act in opposition to each other, q along (1, 1), and b along (−1, 1). We take the same initial condition as in previous numerical experiments.
From the time series in Figure 6, we can see two regimes. The first regime is where the anisotropic diffusion dominates, and we see the orientation of the cells along the positive diagonal. The second regime occurs at the end of the simulation, where the higher density regions are pulled towards the (−1, 1) diagonal, and the lower density regions continue on their path. The split in dominance occurs because the chemotaxis terms' strength is dependent on the density and time. Since the chemoattractant concentration starts at zero, the anisotropic diffusion initially takes control. This state of affairs continues until the anisotropic diffusion squeezes the solution enough to reach a cell density critical point. The chemotaxis terms dominate for intermediate cell densities inducing movement along the orthogonal diagonal. The result is this crisscrossing knitting pattern. To investigate this effect further, we look at the magnitude of the chemotactic and anisotropic fluxes, given by The chemotactic flux F S is proportional to ρ(1 − ρ); hence, it is maximal for ρ = 1/2. It is also proportional to the chemotactic gradient ∇S. The anisotropic flux F q has no ρ dependence in the amplitude, but a ρ dependence inside the divergence term. We plot |F q (x, t)| and |F S (x, t)| for three time points in Figure 7 for the same parameter values as in Figure 6. Initially (at t = 0), we see the initial condition for the anisotropy q in the first row, and since there is no chemical concentration S yet, F S = 0. Note that the first row in Figure 7 only shows the magnitude of F q . The dominant direction is always the diagonal (1, 1). At the time t = 0.0262, we are in the regime where F q continues to diffuse in mainly the (1, 1) direction, whereas F S has just appeared. We see already torus-like structures, which indicate that the chemotactic flux F S is maximal for an intermediate cell density 1/2 and declines again as ρ ≈ 1. As we continue to time t = 0.32, the amplitude of the two fluxes are essentially the same. However, the anisotropic flux is still in the dominant direction (1, 1), while the chemotactic flux is in direction (−1, 1). Hence, we see the competition in those two directions. We see numerically fascinating patterns that have, to our knowledge, never been reported in a chemotaxis system. The effects of anisotropy and chemotaxis interact in complicated ways.

Entrapment by Fibres
For this numerical simulation, we are interested if we can trap a group of cells behind a fiber wall, with some food on the other side to incentivize attempts to escape. We take a horizontally aligned fiber wall as q(v) = 1 4π I 0 (k(x, y)) e k(x,y)v·u + e −k(x,y)v·u , b = 1/50, with u = (1, 0).
The variance matrices (35) and (40) take the form We consider a domain of [0, 250] 2 and assume the linear feature is concentrated at y = 125 with concentration parameter k(x, y) = 5e This choice of k(x, y) causes anisotropy on the line centered at y = 125 along the x-axis and isotropic in other directions. For initial conditions, we take the same Gaussian bumps as before but only place them above the fiber (Figure 8a), and we place a large concentration of the chemoattractant below the fiber wall (Figure 8b). We want to see the cells can transfer through this barrier and find the chemotactic signal. From Figure 8c, we can see that most of the cells continue their normal development as seen in the other simulations, and only a small portion of the cells gets attracted to the fiber wall. They start moving along the fiber wall while a fraction (3.76%) manages to cross and aggregate below the fiber wall. If we increase the amount of chemoattractant, we see more escape, and if we decrease the amount, fewer cells cross the fiber wall. In Figure 8d, we included chemotactic mixing of gradient and anisotropy. For the chemotactic mixing case, we took (18) with the same matrix A, and with s = 3/100π. In that case, the cell movement is aligned in the diagonal direction, and the crossing appears more effective, as 4.3% of the population can cross.

Discussion
In this paper, we analyzed two effects that can lead to directional orientation of moving particles: chemotaxis and anisotropy. The influence of chemotaxis on particle movement, particle aggregations, pattern formation and blow-up has been studied extensively over the past decades. Comprehensive reviews can be found in Horstmann 2003 [48,49], Hillen and Painter 2009 [14], and Bellomo et al. 2015 [16]. If chemotaxis is combined with tissue anisotropy, then new effects arise. For example, Kumar et al. [19] show that the combination of repulsive Ph-taxis and tissue anisotropy can lead to palisade formation, as found in glioma histologies. We show here that anisotropic chemotaxis can explain vasculogenesis as in experiments of Serini et al. [6], that it produces new types of criss-cross patterns, and that it can be used to study the effect of trapping of chemotactic populations through linear features of the tissue. The anisotropic chemotaxis model (34) was derived from a mesoscopic description using kinetic transport equations. This framework includes cell-specific characteristics, such as cell velocity, turning rates, tissue anisotropy, and chemotactic mixing. All these terms are informed from individual particle tracking. Finally, through a parabolic scaling technique, the macroscopic anisotropic chemotaxis model (34) is derived. The numerical solution of the anisotropic chemotaxis model shows new challenges through the presence of the anisotropic terms. To address this challenge, we use a powerful reaction-advectiondiffusion solver of Chertock et al. [21] and modify it to include anisotropy. We can do this without the inclusion of additional diagonal grid points (as done in Kumar [19], for example), rather, using a smart averaging over grid points that are already defined (see the details in Section 5). Using this numerical scheme, we went through numerous simulations; the first numerical experiment we did was to recreate the vascular network assembly results in Serini et al. [6]. We were able to match our solutions to theirs qualitatively; our ability to do this suggests that the flux rapidly relaxes to the equilibrium in their hyperbolic model (2). For the other numerical experiments, we explored the effects of anisotropic diffusion and chemotactic velocity mixing. For the anisotropic diffusion, the effect is essentially smearing the cells along the fiber direction since movement along those directions is favored. The chemotactic velocity had the opposite effect. The mixing dissuaded movement along those directions, but unaffected diffusion along its complement, retrieving similar results as the anisotropic diffusion. Subsequently, we experimented numerically on the inclusion of both the anisotropic diffusion and the chemotactic mixing, which lead to a mixed regime setup. In regions with low cell density, the anisotropic diffusion would dominate, but the chemotactic velocity would govern the motion in high-density regions. Finally, we were interested if chemotactic mixing would help cells move past a fiber wall. We found that the population gets partially trapped by linear tissue features.
The following steps ask for a complete analysis of the observed pattern formation mechanism. It is unknown what type of bifurcation would lead to network patterns and/or criss-cross patterns. It is our plan for future research to investigate the underlying instabilities.
Overall, we show that tissue anisotropy can profoundly affect chemotactic orientation, explaining new biological phenomena.