Direction-Dependent Turning Leads to Anisotropic Diffusion and Persistence
DDirection-Dependent Turning Leads toAnisotropic Diffusion and Persistence
Nadia Loy, ∗ Thomas Hillen, † Kevin J. Painter ‡ January 13, 2021
Abstract
Cells and organisms follow aligned structures in their environment, a process that cangenerate persistent migration paths. Kinetic transport equations are a popular modellingtool for describing biological movements at the mesoscopic level, yet their formulationsusually assume a constant turning rate. Here we relax this simplification, extending toinclude a turning rate that varies according to the anisotropy of a heterogeneous envi-ronment. We extend known methods of parabolic and hyperbolic scaling and apply theresults to cell movement on micro-patterned domains. We show that inclusion of orien-tation dependence in the turning rate can lead to persistence of motion in an otherwisefully symmetric environment, and generate enhanced diffusion in structured domains.
Keywords
Cell migration, Boltzmann equation, persistence, direction dependent turningrate, macroscopic limits.
Subject class (MSC 2020)
Movement of cells through tissues is critical during both healthy and pathological processes.Embryonic development relies on cells migrating from origin to final tissue destination, repairprocesses necessitate movement of fibroblasts and macrophages into the wound site, and mi-gration of cancerous cells, unhappily, leads to tumour invasion and metastasis dissemination.Consequently, there is clear reason to understand the factors that guide cells with one suchprocess, contact guidance , defining the movement of cells along linear/aligned tissue features,for example blood vessels, white matter brain fibres, or the collagen fibres of connective tissue.The influence of contact guidance on cell movement has been considered via a variety ofmathematical approaches [13, 12, 44], with kinetic transport equations proving particularlypopular [22, 34, 24, 6]. Transport equations account for the microscopic features of movement,describing a migration path according to its statistical properties (turning rate, movementspeed and movement direction), and such models for contact guidance have been successfullyapplied to, for example, glioma invasion [35, 26, 48]. Yet these studies have been simplified ∗ DISMA, Politecnico di Torino, Torino, Italy, [email protected] † University of Alberta, [email protected] ‡ DIST, Politecnico di Torino, Torino, Italy, [email protected] a r X i v : . [ q - b i o . CB ] J a n hrough taking turning rates to be independent of orientation, whereas experiments indicateconsiderably more complexity (e.g. [42, 41, 40].Here we extend the known theory of scaling limits for transport equations in biology tocases where the turning rate is direction-dependent. In the physical context, our model canbe seen as a non-homogeneous linear Boltzmann equation with a micro-reversible process inwhich the cross-section is factorised into a turning kernel and a direction-dependent turningrate [38, 11]. The direction dependence of the turning rate augmentation presents mathemat-ical challenges, requiring reflection on how the involved function spaces should be modifiedfor a Fredholm alternative argument to be constructed. We obtain expressions for the macro-scopic diffusion and advection that incorporates the influence of sophisticated turning ratechoices. We apply the model to the movement data of cells on oriented microfabricated sur-faces generated by Doyle et al. [14], showing that strong alignment can lead to persistence ofmovement and, at a macroscopic level, enhanced diffusion.The outline of this paper is as follows. We use the remainder of this introduction toprovide background on contact guidance and detail experimental investigations of cell move-ments on micro-pattern domains. We also review pertinent modelling literature, particularlyusing kinetic transport equations. In Section 2, we introduce the model, explain the basicassumptions, and introduce statistical meaningful quantities such as mean velocity, runtimesalong fibres, directional variance, and persistence. In Section 3 we consider two scaling lim-its, the parabolic limit (Theorem 3.1) and the hyperbolic limit . In particular, we generalisethe technique proposed in [22] to the case of a direction-dependent turning rate and ob-tain a macroscopic diffusive limit with a distinctive structure. Section 4 is used to discusspertinent special cases, with Section 5 illustrating how the new dynamics can result froma direction-dependent turning rate. Extending the analysis to a particularly relevant form,Section 6 is used to demonstrate the utility of the model for describing cell migration pathson microfabricated surfaces. We close with a discussion in Section 7. The extracellular matrix (ECM) is a fundamental ingredient of connective tissues and consti-tutes the major non-cellular component of tissues and organs. Cell migration through ECMcan occur individually or collectively, with individual further classified into amoeboid andmesenchymal forms [50]. Mesenchymal migration is typically slower, with a cell secretingdegrading enzymes (e.g. MMPs) that create space for movement. Thus, mesenchymal migra-tion can significantly alter the local ECM structure. Amoeboid migration is often faster, withfrequent turns and shape changes allowing a cell to squeeze through matrix gaps; contactsare fleeting, leading to moderate and transient changes to the ECM architecture. Cells mayswitch between migration modes, for example in response to the biomechanical resistanceof the ECM. This mesenchymal-amoeboid transition [53] may potentially optimise tumourinvasion [21] in complex heterogeneous micro-environments [49].Regardless of migration type, the architecture of the ECM is a major determinant of move-ment. ECM is formed from various proteins, with collagen often the principal constituent[2]. Individual collagen proteins are organised into cable-like fibres, collectively creating anetwork. Adhesive attachments between cells and matrix-binding sites anchor the cell andprovide the focal points for exerting the forces needed for forward propulsion. Consequently,by protruding and pulling itself along fibres, cells follow the local topology of the matrix( contact guidance ) [15]. The mesh formed from collagen therefore offers an example of a bidi- rectional anisotropic network, bidirectional in the sense that movement preferentially followsfibres but no specific direction is favoured. Anisotropic bidirectional tissues extend to otherenvironments, a particular relevant example being the brain’s white matter. Here it is thelong and bundled neuronal axons that generate the network and its arrangement is believedto be a key determinant in the anisotropic invasion of gliomas [18, 19].The question of how an anisotropic environment influences cell migration is highly suitedto modelling and various approaches have been developed. Agent-based models that incor-porate contact guidance include lattice-free particle approaches (e.g. [10, 30, 43]) and thosebased on the Cellular Potts Model (e.g. [46, 44, 45]) and other automata (e.g. [49]); theindividual-level description is clearly advantageous for incorporating microscopic structure.Continuous models, though, have also been developed, for example the anisotropic biphasictheory (ABT) developed in [13] and transport equations studied in [12]. A transport equa-tion developed in [22] describes the contact-guided migration of mesenchymal (and amoeboid)cells in evolving anisotropic networks of unidirectional or bidirectional type, with this modelextended and subjected to numerical exploration in [34]. Transport models have a ‘stepping-stone’ nature, lying at a point between an individual and macroscopic model: they sit at amesoscopic level, describing the statistical distribution of the individual microscopic velocitiesand positions through density distribution functions. Subsequent up-scaling can generate afully macroscopic model, typically of drift-diffusion nature, and capable of capturing move-ment at a large-tissue level. In the study of [22] the author employed such scaling techniquesto recover macroscopic limits.The transport equation in [22] is predicated on an underlying stochastic velocity-jumpmodel of migration [32], i.e. fixed-velocity runs interspersed with velocity changes. Thetransitional probability for switching velocities (from a pre-reorientation to post-reorientationvelocity) can be decomposed into two elements: a turning rate function that dictates the3ate at which switches occur, and a turning kernel that describes the selection of the newvelocity/direction. The latter was taken to be an angular distribution (potentially spaceand time varying) that encodes the oriented ECM network structure. Thus, contact-guidedmigration was included through an increased likelihood of a cell choosing the dominating localfibre orientation. Consequently, the model captures the anisotropic spread of a population inan aligned bidirectional network and simulations in [34] demonstrate that the environment cansubstantially impact on spatial structuring, for example trapping populations inside or outsideregions of high anisotropy or dictating pathways of invasion. Various real world applicationshave been considered, including predicting the spatial spread of glioma (e.g. [35]) or wolfmovement along seismic lines (e.g. [26]).The turning rate in [22], however, was taken to be independent of orientation. To un-derstand a consequence of this simplification, consider the migration paths of cells subjectedto manufactured environments, such as surfaces subjected to micropatterning (e.g. [14, 51],see Figure 1) or constructed anisotropic collagen networks [42, 41]. In [14] the precise mi-cropatterning of fibronectin on a two-dimensional surface enabled fabrication of controlledanisotropic environments, with the schematics in Figure 1 (a–b) demonstrating two such ar-rangements. Here, cells can move from an effectively 1D region (stripes, replicating highlyaligned parallel fibres) to a 2D region where the 2D regions are both isotropic, but either (a) uniformly isotropic , or (b) featuring criss-crossing perpendicular stripes. As we will explicitlyshow in Section 6, the earlier transport model of [22] is unable to discriminate between thesescenarios.If, instead, cells modulate their turning frequency by turning infrequently when movingalong fibres, we can expect very distinct behaviours. Under the uniformly isotropic case, cellswould be expected to meander significantly in the uniformly isotropic region, adopting shortruns in any orientation. On the other hand, criss-crossing stripes could allow significant trans-lations in either of the two dominating axial directions, hastening rediscovery of the quasi-1Dregions. Indeed, evidence is found of this in the experiments of [14], cf.
Supplementary Movie7. We will explicitly show in Section 6 that direction-dependent turning can lead to enhancedeffective diffusion.
Let p = p ( t, x, w ) denote the cell density distribution, defined at time t ≥
0, position x ∈ Ω ⊆ R d and velocity w ∈ V . We typically assume V to be a compact set as in [24], and inparticular consider V = [ s , s ] × S d − ([22]), where S d − is the set of all possible directionsˆ w ∈ R d and the hat-symbol indicates unit vectors. The limits s , s denote the minimaland maximal speed ( | w | ) of the cells, with 0 ≤ s ≤ s < ∞ . Note that if the speed isapproximately constant we simply set V = s S d − .The governing transport equation for describing cell movement is ∂p ( t, x, w ) ∂t + w · ∇ p ( t, x, w ) = L p ( t, x, w ) , (1)where the operator ∇ denotes the spatial gradient. The turning operator L is a linear operatorthat models the change in velocity of individuals per unit of time at ( x, w ) that is not due to4he free particle drift. L is generally defined as an integral operator on L spaces [24]: L : L ( V ) (cid:55)−→ L ( V ) ,p ( t, x, w ) (cid:55)−→ L p ( t, x, w ) , where ( t, x ) are independent parameters and L p ( t, x, w ) = − µ ( t, x, w ) p ( t, x, w ) + (cid:90) V µ ( t, x, w (cid:48) ) q ( t, x, w, w (cid:48) ) p ( t, x, w (cid:48) ) dw (cid:48) . (2)This operator describes the velocity scattering. As noted earlier, the key determinants of themigration path are the turning rate function, µ ( t, x, w ), and the turning kernel , q ( t, x, w, w (cid:48) ).Viewed in this light, the first term on the right hand side of (2) models particles switchingaway from velocity w and the second one takes into account the particles switching intovelocity w from all other velocities.The turning kernel, q ( t, x, w, w (cid:48) ), denotes the probability measure of switching velocityfrom w (cid:48) to w , given that a turn occurs at location x and time t . Here we adopt the samereasoning as in [22], by assuming reorientation is dominated by the fibrous/anisotropic envi-ronmental structure, such as collagen matrix fibres or white matter tracts. Then, the choiceof new direction is derived from this structure, rather than the incoming velocity, and for q we assume: A1 q ( t, x, w, w (cid:48) ) = q ( t, x, w ) depends only on the post-turning orientation; A2 q ( t, x, w ) ≥ ∀ w ∈ V, a.e. x ∈ Ω , ∀ t ≥ A3 q ( t, x, · ) ∈ L ( V ) and (cid:90) V q ( t, x, w ) dw = 1 a.e. x ∈ Ω , ∀ t ≥ q ( t, x, ˆ w ), defined on R + × Ω × S d − and satisfying˜ q ≥ , (cid:90) S d − ˜ qd ˆ w = 1, to the turning kernel q ( t, x, w ) = ˜ q ( t, x, ˆ w ) ω , ω = (cid:90) V ˜ q ( t, x, ˆ w ) dw. Note that q ( t, x, w ) assumes w ∈ V while ˜ q ( t, x, ˆ w ) is defined for unit vectors only, but theironly difference lies in a constant scaling factor that accounts for the difference between V and S d − . Given this, we will interchangeably call q the turning kernel or the distribution of fibreorientations. As a consequence of A1 - A3 , the operator (2) simplifies to L p ( t, x, w ) = − µ ( t, x, w ) p ( t, x, w ) + q ( t, x, w ) (cid:90) V µ ( t, x, w (cid:48) ) p ( t, x, w (cid:48) ) dw (cid:48) . (3)The turning rate function, µ ( t, x, w ), gives the rate at which velocity switches are madefor a particle located at x at time t , moving in direction w . It is at this point where wesubstantially diverge from [22], lifting the assumptions on µ and allowing w -, x -, and t -dependence. Significantly, this allows the turning rate µ to depend directly on the fibreorientation q , for example allowing a cell to continue movement with the same velocity if it ismoving in the direction of highly aligned fibres. Note that as µ = µ ( t, x, w ) is a turning rate,1 /µ ( t, x, w ) defines the mean time spent by a cell running along a linear tract with velocity w between two consecutive turns performed at time t , location x . We assume:5 µ ( t, x, · ) ∈ L ( V ) , ∀ t > , x ∈ Ω ; M2 q ( t, x, · ) µ ( t, x, . ) ∈ L ( V ) , ∀ t > , x ∈ Ω . To analyse (1) under (3) we make use of a number of statistical properties of the correspondingfibre and turning distributions, such as expectations and variances. A summary of theseexpressions is given in Table 1.1.
Distribution of new directions.
We consider the distribution of newly chosen direc-tions, q , with expectation E q ( t, x ) = (cid:90) V q ( t, x, w ) w dw. (4)This is also the the mean new velocity after a turn and has the variance-covariancematrix V q ( t, x ) = (cid:90) V q ( t, x, w ) ( w − E q ) ⊗ ( w − E q ) dw . (5)2. Cell mean velocity and variance.
We introduce similar macroscopic quantities forthe cell population, although we stress p is not itself a probability measure. First wedefine the macroscopic density of the population p at time t and position x as¯ p ( t, x ) = (cid:90) V p ( t, x, w ) dw , (6)and the total mass of the population in Ω, m ( t ) = (cid:90) Ω ¯ p ( t, x ) dx . Note that, with no population kinetics and assuming suitably lossless boundary con-ditions, the total mass will be conserved in time. With these definitions in place wecan introduce the moments of the normalized cell distribution ˆ p ( t, x, w ) = p ( t, x, w )¯ p ( t, x ) ,which will be a probability measure for all t and x . In particular we can introduce theexpectation E ˆ p ( t, x ) = (cid:90) V ˆ p ( t, x, w ) wdw , which is the mean velocity of the normalized population, and the variance (variance-covariance matrix) V ˆ p ( t, x ) = (cid:90) V ˆ p ( t, x, w ) ( w − E ˆ p ) ⊗ ( w − E ˆ p ) dw . The latter provides information on the width of the distribution ˆ p in different directions.This tensor is symmetric, but can be anisotropic, i.e. the level sets of ˆ w (cid:55)→ ˆ w T V ˆ p ˆ w areellipsoids. With this we can identify the mean velocity of the cell population as (cid:90) V p ( t, x, w ) wdw = ¯ p ( t, x ) E ˆ p ( t, x )6nd the variance-covariance of the population velocity as (cid:90) V p ( t, x, w ) ( w − E p ) ⊗ ( w − E p ) dw = ¯ p ( t, x ) V ˆ p ( t, x ) . Turning part of the population.
The turning operator definition (3) reveals a newmacroscopic quantity, p µ ( t, x, w ) = µ ( t, x, w ) p ( t, x, w ) . (7) p µ can be interpreted as the part of the cell population that moves in direction w andis currently turning . Then, the total turning population per unit time is¯ p µ ( t, x ) = (cid:90) V µ ( t, x, w ) p ( t, x, w ) dw, (8)which is the expression in (3). The turning operator (3) can then be re-written as L p ( t, x, w ) = ¯ p µ ( t, x ) q ( x, w ) − µ ( t, x, w ) p ( t, x, w ) . (9)By normalising p µ we can define the mean incoming velocity of the turning population as E p µ ( t, x ) = (cid:90) V µ ( t, x, w ) p ( t, x, w ) wdw ¯ p µ ( t, x ) , and its variance-covariance matrix accordingly.4. Run times along the fibres.
We discover later that the stationary distributions areproportional to the ratio q/µ . In fact, the corresponding distribution arises in manyforthcoming calculations. Hence, we introduce C ( t, x ) = (cid:90) V q ( t, x, w ) µ ( t, x, w ) dw and the normalised distribution T ( t, x, w ) = 1 C ( t, x ) q ( t, x, w ) µ ( t, x, w ) . (10)Since µ ( t, x, w ) is a turning rate, the function τ ( t, x, w ) = 1 µ ( t, x, w )is the mean time spent moving in direction w . Then, T ( t, x, w ) is the distribution ofrun times along the fibres of the network and its expectation, E T ( t, x ) = 1 C ( t, x ) (cid:90) V q ( t, x, w ) µ ( t, x, w ) w dw , (11)is the average velocity along the fibre distribution. Note that this quantity discriminatesbetween the choice of a parabolic scaling (leading to a diffusive limit, for E T ≈
0) or a7yperbolic scaling (for E T (cid:54) = 0). With this interpretation we can regard the followingnormalisation constant C ( t, x ) = (cid:90) V τ ( t, x, w ) q ( t, x, w ) dw as the mean run time between turns. We can further define the displacement vectorand the mean displacement vector, respectively χ ( t, x, w ) = wτ ( t, x, w ) , ¯ χ ( t, x ) = (cid:90) V χ ( t, x, w ) q ( t, x, w ) dw , (12)such that E T ( t, x ) = ¯ χ ( t, x ) C ( t, x )becomes the ratio of the mean displacement vector over the mean run time on the fibrenetwork.5. Persistence.
Persistence , ψ d , is a measure of a random walker’s tendency to maintaindirection during directional changes. It has values ψ d ∈ [ − , ψ d = 1 denotesperfect persistence (continuing with the previous direction), ψ d = 0 denotes uniformturning and ψ d = − mean cosine along a particle trajectory, however in our abstractframework we define it as the mean velocity of the equilibrium distribution T , i.e. ψ d ( w ) = E T · w | E T || w | . (13)Explicit use of the vector product shows that ψ d does indeed arise as a mean cosine, as ψ d = cos( (cid:94) ( E T , w )) , where (cid:94) ( E T , w ) denotes the angle between E T and w . We have not yet shown that T is the equilibrium distribution, but do so in the next Section. Observe that if µ doesnot depend on w , Eq. (13) becomes ψ d ( w ) = E q · w | E q || w | , (14)as defined in [32]. Hence, for turning rates that do not depend on the direction, thepersistence will vanish for a bi-directional tissue ( i . e . E q = 0). Our extension to w -dependence in µ lifts this limitation, allowing non-vanishing persistence even for a bi-directional tissue with E T (cid:54) = 0. In other words, Eq. (1) with (3) is capable of generatinga persistent random walk even under a fully symmetric configuration.We return to the turning operator (9). As expected, due to assumption A2 , we observe thatthe total cell density during turning will be conserved: (cid:90) V L p ( t, x, w ) dw = (cid:90) V ¯ p µ ( t, x ) q ( t, x, w ) dw − (cid:90) V µ ( t, x, w ) p ( t, x, w ) dw = ¯ p µ ( t, x ) − ¯ p µ ( t, x )= 0 . . distribution meaning expectation variance˜ q ( t, x, ˆ w ) distribution of fibre orientations q ( t, x, w ) distribution of newly chosen directions E q ( t, x ) V q ( t, x )¯ p ( t, x ) = (cid:82) p ( t, x, w ) dw macroscopic cell densityˆ p ( t, x, w ) = p ( t,x,w )¯ p ( t,x ) normalized cell density E ˆ p ( t, x ) V ˆ p ( t, x ) p µ ( t, x, w ) = µ ( t, x, w ) p ( t, x, w ) turning part of the population E p µ ( t, x ) V p µ ( t, x ) T ( t, x, w ) = C ( t,x ) q ( t,x,w ) µ ( t,x,w ) distribution of run times along the fibres E T ( t, x ) V T ( t, x ) The average outgoing velocity of the total turning population, meanwhile, changes: (cid:90) V L p ( t, x, w ) w dw = (cid:90) V ¯ p µ ( t, x ) q ( t, x, w ) w dw − (cid:90) V µ ( t, x, w ) p ( t, x, w ) w dw = ¯ p µ ( t, x ) E q ( t, x ) − ¯ p µ ( t, x ) E p µ ( t, x )= ¯ p µ ( t, x ) (cid:0) E q ( t, x ) − E p µ ( t, x ) (cid:1) . (15)Therefore, the mean velocity of the turning population, E p µ , relaxes towards the averagepost-turning velocity, E q , imposed by the fibre network.For much of the analysis we assume fibre orientations and turning rates are time indepen-dent, i.e. q ( x, w ) and µ ( x, w ). Similar arguments apply for the time dependent case, yet thenotation becomes clumsier. As a first step towards finding equilibrium distributions, we compute the kernel of the turingoperator L . A function φ ( w ) belongs to ker( L ) if and only if it satisfies L φ ( w ) = 0 ⇐⇒ − µ ( x, w ) φ ( w ) + q ( x, w ) (cid:90) V µ ( x, w (cid:48) ) φ ( w (cid:48) ) dw (cid:48) = 0 ⇐⇒ µ ( x, w ) φ ( w ) = q ( x, w ) (cid:90) V µ ( x, w (cid:48) ) φ ( w (cid:48) ) dw (cid:48) . We can write φ ( w ) = q ( x, w ) µ ( x, w (cid:48) ) ¯ φ µ ( x ) , where ¯ φ µ is a w -independent function¯ φ µ ( x ) = (cid:90) V µ ( x, w ) φ ( w ) dw, which is the fraction of turning cells of the population ¯ φ . The w -dependence of elements ofker( L ) is given by the ratio q ( x, w ) /µ ( x, w ), i.e. it is given by the run-time distribution alongthe fibres, T ( x, w ), which we introduced in (10). Then φ ( x, w ) = ¯ φ ( x ) T ( x, w ) , and ker( L ) = (cid:104) T ( x, w ) (cid:105) . Hence a stationary state of the equation (1)-(3) has the form M ( x, w ) = β ( x ) T ( x, w ) (16)9hich we call the Maxwellian. We remark that (1) with (3) is the linear Boltzmann equation.In particular, the quantity σ ( t, x, w, w (cid:48) ) = q ( t, x, w ) µ ( t, x, w (cid:48) )is the cross section and the equilibrium probability density (normalized to 1) given by L ( T ) =0 is (10). The function T will be nonnegative, since q and µ are nonnegative, and it is an L ( V ) function via assumption M2 . Therefore, the stochastic process ruled by σ satisfiesmicro-reversibility, i.e. σ ( x, w, w (cid:48) ) T ( x, w (cid:48) ) = σ ( x, w (cid:48) , w ) T ( x, w ) . Consequently, existence and uniqueness of a nonnegative solution p ∈ L (Ω × V ) of (1),(3)with initial condition p ∈ L (Ω × V ) and non-absorbing boundary conditions [5] is a classicalresult of kinetic theory (see, for example, [37]). We in particular consider a special class ofnon-absorbing boundary conditions, given by no-flux boundary conditions [39].Arguing as in [38, 3], we can prove a linear version of the classical H-Theorem for thelinear Boltzmann equation (1), (3) with initial condition p = p (0 , x, w ) ∈ L (Ω × V ). Let usintroduce the entropy S = − H where, for any given convex function Φ : R + → R + , H Φ [ p | M ]( t ) = (cid:90) Ω (cid:90) V p ( t, x, w )Φ (cid:18) p ( t, x, w ) M ( x, w ) (cid:19) dw dx, p ( t, · , · ) ∈ L (Ω × V ) . In the above, M is the Maxwellian satisfying (cid:90) Ω × V M ( x, w ) dxdw = (cid:90) Ω ¯ p ( x ) dx where we assume that p has finite mass (cid:90) Ω ¯ p dx and entropy H Φ [ p | M ](0). Hence, β ( x ) issuch that (cid:90) Ω β ( x ) dx = (cid:90) Ω ¯ p ( x ) dx and, in particular, if a stationary state p ∞ ( x, w ) exists,then β ( x ) = ¯ p ∞ ( x ), so that p ∞ ( x, w ) = ¯ p ∞ ( x ) T ( x, w ) . Under non-absorbing boundary conditions, via exactly the same procedure followed in [38] itis possible to prove that dH Φ [ p | M ] dt ( t ) ≤ t ≥ dH Φ [ p | M ] dt ( t ) = 0 , iff p ( t, x, w ) = M ( x, w ) , where p ( t, x, w ) is the unique L solution to (1),(3) with initial condition p and non-absorbingboundary conditions. Furthermore, continuing to argue as in [38], it can be proved thatprovided (cid:90) Ω (cid:90) V (cid:0) w + | log p | (cid:1) dwdx < ∞ , thenlim t →∞ (cid:90) Ω (cid:90) V | p ( t, x, w ) − M ( x, w ) | dw dx = 0 . Macroscopic limits
In this section we consider two distinct scalings for the transport equation: i ) the parabolicscaling , and ii ) the hyperbolic scaling . The former applies in a diffusion-dominated case whilethe latter corresponds to the drift-dominated case. The cases differ through the relativescaling of time and space in a suitably small parameter ε .In the parabolic scaling we consider a small parameter, ε (cid:28)
1, and assume macroscopictime and space scales ( τ, X ) that scale according to τ = ε t, X = εx. A paradigm for scaling in this manner can be found by comparing the microscopic scales of runand tumble movements of
E. coli bacteria and the experimental scales at which population-level phenomena form, such as travelling bands [1] or cellular aggregates [4]. Runs are charac-terised with a movement speed typically around 10-20 µm/s with a tumble taking place everysecond or so. Large scale patterning phenomena typically arise after a few hours, i.e. O (10 )seconds. Hence, ε = mean run time/experimental time = 10 − and, in turn, ε = 10 − .Given the micron spatial scale of individual movement, the corresponding macroscopic scaleis 10 × µm = 1 mm , which is the scale of the cell aggregates that typically form.This scaling therefore demands a sufficiently slow time scale over which diffusion can beginto dominate or, equivalently, when cells have large speeds and turning rates. Rescaling (1),(3)gives ε ∂p∂τ ( τ, X, w ) + εw · ∇ p ( τ, X, w ) = − µ ( X, w ) p ( τ, X, w ) + q ( X, w ) (cid:90) V µ ( X, w (cid:48) ) p ( τ, X, w (cid:48) ) dw (cid:48) , (17)where the gradient ∇ is now applied with respect to X . We assume that the rescaled turningfrequency and kernel are such that µ ( X, w ) , q ( X, w ) ∼ O (1). Thanks to classical results (seee.g. [37]), we have existence and uniqueness of the solution of (17) with initial condition p ∈ L (Ω × V ) and non-absorbing boundary conditions. We now consider an asymptoticexpansion of p in orders of ε , p ( τ, X, w ) = p ( τ, X, w ) + εp ( τ, X, w ) + ε p ( τ, X, w ) + O ( ε ) , (18)and we are particularly interested in the leading order term p . In the diffusion-dominated case we assume that the macroscopic drift term E T = 0. Weconsider this case first to introduce the necessary technical notations. Theorem 3.1.
Let assumptions A1 - A3 be satisfied and assume E T = 0 . (19) Consider equation (17) with expansion (18). The leading order term, p ( τ, X, w ) , of equation(18) satisfies p ( τ, X, w ) = ¯ p ( τ, X ) T ( X, w ) , here ¯ p ( τ, X ) solves the macroscopic anisotropic diffusion equation ∂∂τ ¯ p ( τ, X ) = ∇ · (cid:90) V µ ( X, w ) ∇ · (cid:104) ¯ p ( τ, X ) D T ( X, w ) (cid:105) dw , (20) with microscopic anisotropic diffusion tensor D T ( X, w ) = w ⊗ w T ( X, w ) . Proof.
We have seen that ker( L ) = (cid:104) T (cid:105) . In order to invert L on the orthogonal complement of its kernel, we use a specific weight forthe inner product in L η ( V ), with η ( X, w ) = µ ( X, w ) T ( X, w ) , (21) i . e . given f, g ∈ L the weighted scalar product is defined by( f, g ) η = (cid:90) V f ( w ) g ( w ) η ( X, w ) dw. We observe that φ ( X, w ) ∈ (cid:104) T (cid:105) ⊥ if and only if0 = (cid:90) V φ ( X, w ) T ( X, w ) µ ( X, w ) T ( X, w ) dw = (cid:90) V φ ( X, w ) µ ( X, w ) dw = ¯ φ µ ( X ) . The range of L is given by all functions that integrate to zero, since (cid:90) V L φ ( X, w ) dw = L φ = 0from mass conservation. Then, this range, as subset of L η ( V ), can be written asRange( L ) = (cid:104) η − (cid:105) ⊥ , since 0 = (cid:90) V φ ( X, w ) dw = (cid:90) V φ ( X, w ) η − ( X, w )( η ( X, w ) dw ) = ( φ, η − ) η . Then L ⊥ : (cid:104) T (cid:105) ⊥ → (cid:104) η − (cid:105) ⊥ , with L ⊥ = L| (cid:104) T (cid:105) ⊥ and this restricted operator L ⊥ is the operator we would like to invert. Given ψ ∈ (cid:104) η − (cid:105) ⊥ , weneed to find φ ∈ (cid:104) T (cid:105) ⊥ such that L φ = ψ . Applying L from (9) we obtain (ignoring argumentsfor clarity of presentation) L φ = ¯ φ µ q − µφ = ψ. Since on (cid:104) T (cid:105) ⊥ we have ¯ φ µ = 0, we can solve for φ as φ ( X, w ) = − µ ( X, w ) ψ ( X, w ) . (cid:16) L ⊥ (cid:17) − : (cid:104) η − (cid:105) ⊥ → (cid:104) T (cid:105) ⊥ , ψ (cid:55)→ − µ ψ. (22)We may observe that the pseudo-inverse depends on the microscopic velocity and on themacroscopic space coordinate through µ ( X, w ).We now substitute expansion (18) into equation (17) and match orders of ε . • For ε we find L p ( τ, X, w ) = 0 , hence p ∈ ker( L ) and we can write p ( τ, X, w ) = ¯ p ( τ, X ) T ( X, w ) . (23) • At order ε we have ∇ · (cid:0) w p ( τ, X, w ) (cid:1) = L p ( τ, X, w ) . (24)To solve this equation for the next order correction term p , we need to invert L . We sawearlier in (22) that L is invertible on (cid:104) η − (cid:105) ⊥ . Hence we check the solvability condition ∇ · ( wp ( τ, X, w )) ∈ (cid:104) η − (cid:105) ⊥ . Using (23) this condition becomes0 = ( ∇ · ( wp ) , η − ) η = (cid:90) V ∇ · ( w ¯ p ( τ, X ) T ( X, w )) η ( X, w ) − η ( X, w ) dw = ∇ · (cid:20)(cid:90) V wT ( X, w ) dw ¯ p ( τ, X ) (cid:21) = ∇ · [ E T ( X )¯ p ( τ, X )] . Hence in this step it is necessary to assume E T ( X ) = 0 . (25)In words, we assume the distribution of run times along the fibre network has no dom-inant direction. In this case we can invert (24) and find p ( τ, X, w ) = 1 µ ( X, w ) ∇ · ( w ¯ p ( τ, X ) T ( X, w )) . (26) • In ε : ∂∂τ p ( τ, X, w ) + ∇ · (cid:0) p ( τ, X, w ) w (cid:1) = L p ( τ, X, w ) . Integrating this equation over V we obtain (cid:90) V ∂∂τ p ( τ, X, w ) dw + ∇ · (cid:90) V wp ( τ, X, w ) dw = 0 . p and (26) for p , we obtain ∂∂τ ¯ p ( τ, X ) (cid:90) V T ( X, w ) dw (cid:124) (cid:123)(cid:122) (cid:125) =1 = ∇ · (cid:90) V µ ( X, w ) w ⊗ w ∇ (cid:104) ¯ p ( τ, X ) T ( X, w ) (cid:105) dw = ∇ · (cid:90) V µ ( X, w ) ∇ · (cid:104) ¯ p ( τ, X ) w ⊗ wT ( X, w ) (cid:105) dw = ∇ · (cid:90) V µ ( X, w ) ∇ · (cid:104) ¯ p ( τ, X ) D T ( X, w ) (cid:105) dw , with an anisotropic diffusion tensor D T ( X, w ) = w ⊗ w T ( X, w ) . D T is a macroscopic scale diffusion tensor, where V T ( X ) = (cid:90) V D T ( X, w ) dw describes the variance of the run time distribution along the network fibres. Lemma 1.
The above parabolic limit equation (20) can be written as an anisotropic drift-diffusion model ∂ ¯ p ∂τ ( τ, X ) + ∇ · (cid:0) a ( X )¯ p ( τ, X ) (cid:1) = ∇ · ∇ · (cid:0) D ( X ) ¯ p ( τ, X ) (cid:1) , (27) with macroscopic diffusion tensor, D , and advection speed, a , given by D ( X ) = (cid:90) V w ⊗ w T ( X, w ) µ ( X, w ) dw (28) a ( X ) = − (cid:90) V w ⊗ w ∇ µ ( X, w ) µ ( X, w ) T ( X, w ) µ ( X, w ) dw. (29) Proof.
The proof relies on straightforward application of the quotient rule. Omitting argu-ments for readability , ∇ · ∇ · (cid:18) ¯ p (cid:90) V w ⊗ w Tµ dw (cid:19) = ∇ · (cid:90) V µ ∇ · (¯ p w ⊗ wT ) dw − ∇ · (cid:18) ¯ p (cid:90) V ∇ µµ w ⊗ wT dw (cid:19) . This representation allows some interesting physicobiological interpretations. • In (21) earlier we defined the weight function η = µT . We can write the above macro-scopic quantities in terms of η as D ( X ) = (cid:90) V w ⊗ w dwη ( X, w ) , a ( X ) = − (cid:90) V w ⊗ w ∇ ln( µ ( X, w )) dwη ( X, w ) . D then appears as an anisotropy matrix related to the measure η − dw , while a measuresthe anisotropic logarithmic gradient with the same measure η − dw .14 We can also relate the terms back to the original network structure given by q ( X, w ).Recall that T = qCµ , where C ( X ) was a normalisation constant. By considering thedistance travelled in direction w we can write D ( X ) = 1 C ( X ) (cid:90) V χ ( X, w ) ⊗ χ ( X, w ) q ( X, w ) dw , a ( X ) = − C ( X ) (cid:90) V χ ( X, w ) ⊗ χ ( X, w ) ∇ ln( µ ( X, w )) q ( X, w ) dw , using the definition of χ from (12). D is then the scaled variance-covariance matrix ofthe directed mean run-time along the fibre network and a is the scaled mean logarithmicderivative of the turning rate, weighted by the mean run times along the fibre network. • The equations simplify drastically when the turning rate µ does not depend on space.In this case ∇ µ = 0 and hence a = 0, i.e., we get a pure anisotropic diffusion equation.Viewed this way, the advection velocity a clearly results from spatial variation in theturning rate, µ ( X, w ), and relates to an anisotropic taxis term measuring the drift dueto the gradient of µ . Spatial dependence in µ ( X, w ) generates an advection pushingcells towards decreasing values of the turning frequency.Due to the fact that V and Ω are compact, we can directly apply a convergence result of[11], giving us Lemma 2.
Suppose A1 - A3 and M1 - M2 hold and assume (19) . Let us also suppose that ∃ C , C ≥ such that | w · ∇ q ( x, w ) | ≤ C µ ( x, w ) , a.e. in Ω × V and (cid:90) V w q ( x, w ) µ ( x, w ) dw ≤ C a.e. in Ω . Let the initial condition p (cid:15) ( x, w ) satisfy (cid:90) Ω × V p (cid:15) ( x, w ) T ( x, w ) dxdw < ∞ and ¯ p (cid:15) (cid:42) ¯ p in H − (Ω) .Let p (cid:15) be the solution to (17) . Then there exists a subsequence ¯ p (cid:15) (cid:42) ¯ p in L ((0 , T ) × Ω) where ¯ p satisfies (27) with initial condition ¯ p ( x ) . The authors of [20] study the parabolic scaling also in the case where the macroscopic drift E T (cid:54) = 0. The key lies in a transformation to moving spatial coordinates, shifting the solutionin the direction of E T as Z = X − E T t . This method also works here and we introduce p ( τ, X, w ) := u ( τ, X − E T τ, w ) . Then, the transport equation (1) transforms to ∂∂τ u + ∇ · [( w − E T ) u ] = L u. Therefore, the scaled transport equation (17) ε ddτ p + ε ∇ · ( wp ) = L p ε ∂∂τ u + ε ∇ · [( w − E T ) u ] = L u. For this modified transport equation we perform the same scaling analysis as before. Weconsider u ( τ, Z, w ) = u ( τ, Z, w ) + εu ( τ, Z, w ) + ε u ( τ, Z, w ) + · · · . Upon comparing orders of ε we find, to leading order, that L u = 0 . Hence, u is in the kernel of L and we can write u ( τ, Z, w ) = ¯ u ( τ, Z ) T ( X, w ) . The order ε terms are ∇ · (cid:2) ( w − E T )¯ u T (cid:3) = L u . To solve this equation for u we require the solvability condition ∇ · (cid:20)(cid:90) V ( w − E T ) T dw ¯ u (cid:21) = 0 , which is true since (cid:90) V ( w − E T ) T dw = E T − E T = 0 . Then u = − µ ∇ · (cid:2) ( w − E T )¯ u T (cid:3) . The terms of order ε are ∂∂τ (¯ u T ) + ∇ · (cid:2) ( w − E T ) u (cid:3) = L u . We integrate this equation over V to obtain0 = ∂∂τ ¯ u − ∇ · (cid:90) V ( w − E T ) (cid:18) µ (cid:19) ∇ · (cid:2) ( w − E T )¯ u T (cid:3) dw (cid:124) (cid:123)(cid:122) (cid:125) ( I ) . (30)Generalising the previous definitions of the diffusion tensor (28) and the drift velocity (29)we define now a macroscopic diffusion tensor D h and macroscopic drift velocity a h as D h ( X ) = (cid:90) V ( w − E T ) ⊗ ( w − E T ) Tµ dw , (31) a h ( X ) = (cid:90) V ( w − E T ) ∇ · (cid:20) w − E T µ (cid:21) T dw . (32)Note that for E T = 0 we return to the previous definitions (28) and (29) of the diffusion-dominated case. 16herefore we have that ∇ · ( D h ¯ u ) = (cid:90) V ( −∇ · E T )( w − E T ) Tµ dw ¯ u + (cid:90) V ( w − E T ) 1 µ ∇ · (cid:2) ( w − E T )¯ u T (cid:3) dw (cid:124) (cid:123)(cid:122) (cid:125) ( I ) + (cid:90) V ( w − E T ) ⊗ ( w − E T ) T ∇ (cid:18) µ (cid:19) dw ¯ u = ( I ) + (cid:90) V ( w − E T ) (cid:20) ∇ · ( w − E T ) 1 µ + ( w − E T ) · ∇ (cid:18) µ (cid:19)(cid:21) T dw ¯ u = ( I ) + a h ¯ u . (33)With these definitions we obtain from (30) a fully anisotropic drift-diffusion model for ¯ u : ∂∂τ ¯ u + ∇ · ( a h ¯ u ) = ∇ · ∇ · ( D h ¯ u ) . (34)Finally, transforming back to ¯ p ( τ, X ) = ¯ u ( τ, X − E T τ ) we find ∂∂τ ¯ p + ∇ · (( a h + E T )¯ p ) = ∇ · ∇ · ( D h ¯ p ) . (35)If E T = 0 we get the same result as already shown in Theorem 3.1 and equation (27) inLemma 1. In drift-dominated phenomena we expect that nondimensionalisation leads to a scaling inwhich the macroscopic time and space scales are of the form τ = εt, X = εx , ε (cid:28)
1. Ef-fectively, cells do not have a large turning frequency, drift dominates and we can perform ahyperbolic limit. The rescaled equation is ε ∂p∂τ + εw · ∇ p = L p . (36)We consider the following expansion p = p + εg + O ( ε ) , (37)where we assume that the correction term g carries zero mass (¯ g = 0) and, therefore, ¯ p = ¯ p .Substituting (37) into (36) we find that the leading order terms ( ε ) are L p = 0, hence p ( τ, X, w ) = ¯ p ( τ, X ) T ( X, w ) . (38)With this choice of p the remaining terms in (36) are ∂p∂τ + εg τ + w · ∇ p + εw · ∇ g = L g + O ( ε ) . (39)We integrate this equation over V , use the form of p from above (38), and assume the O ( ε )-terms are negligible. Then ∂ ¯ p∂τ + ε (cid:90) V g τ dw + ∇ · (cid:90) V w ¯ pT dw + ε ∇ · (cid:90) V wgdw = 0 . g carries no mass, we have ∂∂τ ¯ g = 0. Using the definition of the expectation of T from(11) we obtain ∂ ¯ p∂τ + ∇ · ( E T ¯ p ) + ε ∇ · (cid:90) V wgdw = 0 , (40)which, to leading order, becomes the pure drift model ∂ ¯ p∂τ ( τ, X ) + ∇ · (cid:16) E T ( X )¯ p ( τ, X ) (cid:17) = 0 . (41)The macroscopic drift velocity E T is the expected movement direction based on the av-erage time spent on the fibres. If E T is of order one then (41) is the leading order model.However, for small or zero expectation E T , we can compute the next-order correction term g . Specifically, we use (cid:90) V wg dw = (cid:90) V ( w − E T ) g dw and rewrite equation (40) as ∂ ¯ p∂τ + ∇ · (cid:0) E T ¯ p (cid:1) = − ε ∇ · (cid:104)(cid:90) V ( w − E T ) gdw (cid:105) . (42)From the previous equation (39) we find to leading order that L g = ∂p∂τ + w · ∇ p = ∂ ¯ p∂τ T + w · ∇ (¯ pT ) = −∇ · ( E T ¯ p ) T + w · ∇ (¯ pT ) , (43)where we used the drift model (41) in the last step. To solve for g we need to satisfy thesolvability condition 0 = (cid:90) V −∇ · ( E T ¯ p ) T + w · ∇ (¯ pT ) dw = −∇ · ( E T ¯ p ) + ∇ · (cid:18)(cid:90) V wT dw ¯ p (cid:19) = −∇ · ( E T ¯ p ) + ∇ · ( E T ¯ p ) , (44)which is true since (cid:82) T dw = 1. Hence the right hand side of equation (43) is in the range of L .To solve an equation of the form L g = R with R ∈ Range( L ), we can use the pseudoinverseof L ⊥ and add a term which lies in the kernel of L , i.e. we write g ( τ, X, w ) = − µ ( X, w ) R ( τ, X, w ) + α ( τ, X ) T ( X, w ) , where α ( τ, X ) is independent of w . Here we assumed the correction term has no mass, ¯ g = 0,and hence we choose α ( τ, X ) in such a way that ¯ g = 0. We solve (43) for the correction term, g , to obtain g = − µ (cid:104) −∇ · ( E T ¯ p ) T + w · ∇ (¯ pT ) (cid:105) + αT with α ( τ, X ) = − (cid:90) V µ (cid:104) −∇ · ( E T ¯ p ) T + w · ∇ (¯ pT ) (cid:105) dw. (45)Note that if the turning rate is independent of velocity w (i.e. µ ( X, w ) = µ ( X )), then from(44) we note α = 0 and return to the standard case discussed in the introduction (see [25]).18e write g in the following form, more convenient for later manipulation: g = − µ (cid:104) T ( w − E T ) · ∇ ¯ p + ¯ p ( w · ∇ T ) − ¯ p (( ∇ · E T ) T ) (cid:105) + αT = − µ (cid:104) T ( w − E T ) · ∇ ¯ p + ∇ T · ( w − E T )¯ p + E T · ∇ T ¯ p + (cid:2) ∇ · ( w − E T ) (cid:3) ¯ pT (cid:105) + αT = − µ (cid:104) ∇ · (cid:2) ( w − E T )¯ pT (cid:105) + E T · ∇ T ¯ p (cid:105) + αT. Then, the term in the square brackets of (42) becomes (cid:90) V ( w − E T ) gdw = − (cid:90) V ( w − E T ) 1 µ ∇ (cid:2) ( w − E T )¯ pT (cid:3) dw (cid:124) (cid:123)(cid:122) (cid:125) ( I ) (46) − (cid:90) V µ ( w − E T ) E T · ∇ T dw ¯ p + (cid:90) V ( w − E T ) αT dw (47)= −∇ · ( D h ¯ p ) + a h ¯ p − (cid:90) V ( w − E T ) E T · ∇ Tµ dw ¯ p . (48)In the above we used the relation (33) between the integral ( I ) and the diffusion and driftterms D h and a h , respectively, as well as the fact that (cid:90) V ( w − E T ) αT dw = α (cid:18)(cid:90) V wT dw − E T (cid:90) V T dw (cid:19) = 0 . Combining these calculations with the macroscopic limit equation (42), we obtain the hyper-bolic limit equation with correction term as ∂ ¯ p∂τ + ∇ · (cid:0) ( E T + ε a h )¯ p (cid:1) = ε ∇ · ∇ · (cid:16) D h ¯ p (cid:17) + ε ∇ · (cid:18)(cid:90) V ( w − E T ) E T · ∇ Tµ dw ¯ p (cid:19) . (49)where D h and a h are given by (31) and (32), respectively.A particularly pleasing aspect of the above hyperbolic limit lies in its generalisation ofthe earlier parabolic limit. Specifically, for the case E T = 0 we obtain the same macroscopicquantities: (31) coincides with (28) and (32) coincides with (29). Moreover, for E T = 0,the rather clumsy integral correction term vanishes and we obtain a rescaled version of theparabolic limit equation (27): ∂ ¯ p∂τ + ε ∇ · (cid:0) ( a h )¯ p (cid:1) = ε ∇ · ∇ · (cid:16) D h ¯ p (cid:17) . (50)A further scaling of time with ε then reproduces the parabolic limit (27). For time-varying tissues, q , µ and C will all depend on time. Rescaling (1),(3) with X = εx ,regardless of using time scale τ = ε t or τ = εt , and comparing equal orders of ε allows us toobtain the leading order function of (18), that is p ( τ, X, w ) = ¯ p ( τ, X ) T ( τ, X, w )19here the equilibrium distribution is now also time dependent, T ( τ, X, w ) = q ( τ, X, w ) µ ( τ, X, w ) .The diffusive limit does not change, but the hyperbolic limit has different correction terms.Proceeding as in the previous section, we find that the correction g is a solution to L g = ∂∂τ (¯ pT ) + w · ∇ (cid:16) ¯ pT (cid:17) + O ( ε ) . (51)Let us suppose again that g is of the form g = − R ( τ, X, w ) µ ( τ, X, w ) + α ( τ, X ) T ( τ, X, w )where again α ( τ, X ) does not depend on w . Inverting (51), we find g = − µ ∂ ¯ p∂τ T + ¯ p ∂T∂τ + w · ∇ (cid:16) ¯ pT (cid:17) + α ( τ, X ) T ( τ, X, w )and then g = − µ (cid:104)(cid:16) wT − E T T (cid:17) · ∇ ¯ p + (cid:16) w · ∇ T − ∇ · E T T + ∂T∂τ (cid:17) ¯ p (cid:105) + α ( τ, X ) T ( τ, X, w ) . Since we assumed g carries no mass and have ∂∂τ ¯ T = 0, we obtain the same α as in (45).Concluding, equation (49) becomes in this case ∂ ¯ p∂τ + ∇ · (cid:0) ¯ p ( E T + (cid:15)a h ) (cid:1) = ε ∇ · ∇ · (cid:16) D h ¯ p (cid:17) + ε ∇ · (cid:18)(cid:90) V ( w − E T ) E T · ∇ Tµ dw ¯ p (cid:19) + ε ∇ · (cid:16) ∂∂τ E T ¯ p (cid:17) . (52) To appreciate the power of the analysis, we use this section to consider special cases relevantto various biological scenarios. First, consider the parabolic limit equation (27). Applicationsoften only consider the macroscopic model for ¯ p , and we simplify the notation by using u ( t, x )instead of the cumbersome ¯ p ( τ, X ). We rewrite the parabolic limit (27) in a more standardFickian form ∂u∂t ( t, x ) = ∇ · (cid:2) D ( x ) ∇ u ( t, x ) (cid:3) + ∇ · (cid:2) ( a ( x ) + ∇ · D ( x )) u ( t, x ) (cid:3) . (53)Therefore, the fully anisotropic advection-diffusion process described by the Fokker-Planckequation (27) is revealed as a standard anisotropic Fickian diffusion with an advective com-ponent given by the combination of drift velocity a ( x ) and the gradient of the diffusion tensor ∇ · D ( x ). We consider some special cases and simplify notation by omitting dependencieswhen they are obvious.
1. Suppose µ does not depend on w . The run time distribution, T , then becomes T ( x, w ) = q ( x, w ) C ( x ) µ ( x ) , with C ( x ) = (cid:90) V q ( x, w ) µ ( x ) dw = 1 µ ( x ) . T ( x, w ) = q ( x, w )and the model reduces to previously studied cases (e.g. see [22, 36]). The diffusionmatrix here is given by the variance-covariance matrix of the fibre network distribution, q ( x, w ), divided by the turning rate D ( x ) = (cid:90) V w ⊗ w T ( x, w ) µ ( x ) dw = 1 C ( x ) µ ( x ) (cid:90) V w ⊗ wq ( x, w ) dw = V q ( x ) µ ( x ) . The drift velocity (29), meanwhile, is a ( x ) = − ∇ µ ( x ) µ ( x ) V q ( x ) . (54)In this case, if q is even and then E T vanishes, we obtain a parabolic limit equation (27)as ∂u∂t − ∇ · (cid:20) V q ∇ µµ u (cid:21) = ∇ · ∇ · (cid:18) V q µ u (cid:19) . (55)This can again be written in Fickian form, where we obtain ∂u∂t = ∇ · (cid:18) V q µ ∇ u (cid:19) + ∇ · (cid:20) ∇ · V q µ u (cid:21) . (56)When µ is spatially heterogeneous, the diffusion process (55) is a fully anisotropicprocess with advection given by the taxis term (54). This leads the dynamics towardsdecreasing values of µ . The turning rate µ also scales the diffusivity, with large diffusivityfor small turning rates and vice versa. Hence, frequently turning particles will not showlarge diffusive spread compared to those turning less often. The drift in the direction ofthe negative gradient of µ indicates a drift away from regions of low diffusion towardsregions of high diffusion.2. If we further assume µ is constant and consider q = q ( x, w ), the Fokker-Planck equation(27) leads to ∂u∂t = 1 µ ∇ · ∇ · ( V q u ) . (57)When the oriented habitat is not spatially homogeneous, (57) describes a fully anisotropicprocess where diffusion is again linked to an anisotropic environment (which determinesthe direction of anisotropy) and to the frequency of reorientations (which determinesthe intensity of diffusion).3. Suppose now µ is constant and further assume q = q ( w ), i.e. the network is spatially ho-mogeneous but can still describe an oriented environment. Then the diffusion equation(57) becomes ∂u∂t = 1 µ ∇ · ( V q ∇ u ) . (58)Equation (58) describes diffusion in a spatially homogeneous environment, where anisotropyis due to the presence of a dominant alignment in the environment.21. Suppose both turning rates and turning kernels depend on w but not x , i.e. µ = µ ( w )and q = q ( w ). If we further assume E T = 0, in this case ∇ µ ( X, w ) = 0, the macroscopicdrift from (29 ) satisfies a = 0 and the macroscopic diffusion (28) is D = (cid:90) V w ⊗ w T ( X, w ) µ ( w ) dw. (59)We use this case later to analyse the movement patterns of cells migrating on fabricatedanisotropic surfaces, such as the fibronectin stripe arrangements devised in [14] (see theschematic in Figure 1 and the simulations in Figure 7).5. Finally, suppose a constant fibre distribution, i.e. q = 1 / π, (60)but µ = µ ( w, x, t ), i.e. the turning rate can depend on the orientation of the envi-ronment. If we suppose E T = 0, we are again in the previous case, but the eventualanisotropy will be the direction orthogonal to the dominant direction of µ , as it is theorientation of cells that tend to turn more frequently. Generally we observe an anisotropic diffusion process (28), a consequence of both cells ori-enting with respect to the network alignment via q ( x, w ) and the direction-dependent turningfrequency µ ( x, w ). Specifically, diffusive anisotropy is encapsulated through the second mo-ment of the distribution T ( x, w ) µ ( x, w ) = q ( x, w ) C ( x ) µ ( x, w ) . The nominator/denominator localisations of q and µ naturally suggest that turning directionand turning rate might have opposing impact on macroscopic dynamics.1. For example, suppose that q ( x, w ) ∼ µ ( x, w ) . (61)Here we have T ( x, w ) ∼ µ ( x, w ) and obtain a macroscopic drift term of the form E T ( x ) ∼ (cid:82) V wµ ( x, w ) dw. The macroscopic diffusion tensor (28) in this case will be isotropic D ∼ (cid:90) V w ⊗ wdw = cI. Hence, for (cid:82) wµ ( x, w ) dw (cid:54) = 0, the macroscopic regime will be isotropic and drift-dominated.2. Suppose instead q ( x, w ) ∼ µ ( x, w ) . (62)Here, T ( x, w ) is constant and the macroscopic regime will always be diffusive with E T = 0. The diffusion might be anisotropic, since D ∼ (cid:90) V w ⊗ w q ( x, w ) dw.
22. It is also relevant to consider the case in which q ( x, w ) ∼ µ ( x, w ) , (63)because in many biological cases cells are likely to both choose the fibre direction andturn less frequently when moving in the direction of dominating alignment. In this casewe have both a drift-driven phenomenon, with E T = (cid:82) w µ ( x,w ) dw , and the diffusioncan be anisotropic according to D ∼ (cid:90) V w ⊗ wq ( x, w ) dw. The somewhat curious diffusive terms derived here arise due to the study of active movers :biological particles, such as cells and organisms, generate their own movement energy andare therefore unconstrained by energy and momentum conservation as demanded by classicalphysics. For passive movers (movers simply transported by a fluid environment or some otherstream) the situation is different: in that context, an equation of the form (58) corresponds tothe equation describing diffusion in an anisotropic but homogeneous environment, i.e. withoutspatial variation; equation (57) corresponds to the equation derived in [7], describing particlediffusion in a heterogeneous environment; equation (55) corresponds to the case in which µ depends on the position, and it is then the only case in which the macroscopic equation hasthat particular mixed structure. Within the diffusion theory of physical particles, this mixedstructure arises when describing a so-called thermal effect, i . e . when there is a temperaturegradient [52]. Here, the variable turning rate can be viewed somewhat analogously to tem-perature in a system of physical particles , and the two factors influencing the dynamics arethe spatial heterogeneity and the adaptability to environmental heterogeneity.An example of adaptability in a biological context is provided in [9, 8], addressing starvation-driven diffusion. Here, starvation induces the organism to increase motility and find a betterenvironment, even if it is not known where that may be, i . e . the migration is unbiased. Wefurther note that the three forms of diffusion equation (55),(57),(58) can be derived fromspace-jump processes, where variation in jumping depends on environment assessment of thecurrent location, between locations or the destination [31, 9, 29], with spatially homogeneousjumping times. When both the jumping time and length are spatially heterogeneous, a mixedstructure similar to (27) arises. We present simulations to illustrate key model features, numerically integrating the kinetictransport equation (1),(9) to approximate the density distribution p and in turn the corre-sponding macroscopic density via (6). For computational convenience we restrict to a 2Dspatial setting, a rectangular 2D region Ω = [0 , L x ] × [0 , L y ], and restrict the velocity spaceto V = S , so that w = ˆ w and the speed is set at unit value. In particular, we integrate In a biological context, this ‘temperature’ should not be thought of in its everyday sense, but the macro-scopic temperature that arises from the movement of particles viewed as a mechanical multi-particle system,like in the thermodynamic theory of gases. p µ has to be computed from the p obtained from the current time step.Initially the population macroscopic density, ¯ p ( x ) = ¯ p (0 , x ), is described by a tightly-concentrated Gaussian distribution centred at location ( x , y ), e.g. see Figure 3(a), while cellorientations are uniformly distributed over S . At the boundaries we set diffusive boundaryconditions [27], yielding no-flux boundary conditions at the macroscopic level [39, 28] both inthe parabolic and hyperbolic limit.The remainder of this section is divided into three core tests, designed to illustrate how q and µ alter the macroscopic dynamics. Test 1 explores the extent to which anisotropic orisotropic behaviour emerges with the relationship between q and µ . In Test 2 we demonstratethe tactic effect induced by the spatial gradient of µ . Test 3 extends the analysis to morecomplicated network structures. We specify q and µ using bimodal von Mises distributions, a standard circular distributionwith known analytical forms for the first and second moments (e.g. see [23]). For notationalconvenience we introduce the following short-hand for a bimodal von-Mises distribution, witha given concentration parameter k > θ : bvM k,θ ( w ) = 14 πI ( k ) (cid:16) e kw · θ + e − kw · θ (cid:17) . (64)To further simplify notation we identify a unit vector by its angle, i.e. writing θ = (cos θ, sin θ ) T .Fibres are arranged in two principal patterns, schematised in Fig. 2. We base q and µ on(64), stipulating functions θ q and k q for q and θ µ and k µ for µ . k q = 0 therefore correspondsto unaligned fibres and large k q generates highly aligned fibres along the axis θ q , θ q + π . Testsare devised as follows. • Test 1 . Here we set θ q = µ q = π/ k q = 50 e − . x − +( y − ) and µ as one of µ ∼ constant, ∼ q, ∼ q , ∼ /q. • Test 2 . We set θ q , µ q and k q as in Test 1, but shift k µ to 50 e − . x − . +( y − . ) . • Test 3 . Setting Ω XY = Ω X ∩ Ω Y , where Ω X = { ( x, y ) : 4 ≤ x ≤ } , Ω Y { ( x, y ) : 4 ≤ y ≤ } , we consider q ( x, w )= π on Ω − Ω XY , bvM k q ,θ q on Ω X , bvM k q ,θ ⊥ q on Ω Y , (cid:16) bvM k q ,θ q + bvM k q ,θ ⊥ q (cid:17) on Ω XY ,µ ( x, w )= π on Ω − Ω XY , bvM k µ ,θ µ on Ω X , bvM k µ ,θ ⊥ µ on Ω Y , (cid:16) bvM k µ ,θ µ + bvM k µ ,θ ⊥ µ (cid:17) on Ω XY . (65)Specifically, we will set θ q = π and k q = 50 to form the cross configuration on the rightof Fig. 2, where away from the cross fibres are isotropic and on the horizontal (vertical)arm fibres are aligned in the horizontal (vertical) direction. At the centre, fibres cross.For µ we consider related but subtly distinct forms, allowing distinct turning rate toturning distribution behaviour. 24igure 2: Left: Arrangement and intensity of fibre orientation ( q ) for Tests 1 & 2 . Themaximum of µ is positioned at the red square for Test 1 and the green circle for Test 2.Right: Arrangement and intensity of fibre orientation ( q ) for Test 3 . According to theprecise simulation, we again shift the position of the maximum of µ from the red square tothe green circle. Test 1 was designed to explore the extent to which isotropic or anisotropic behaviour ariseswith the choices of q and µ , for example related as in Eq. (61), (62) or (63). Initialisingaccording to Fig. 3(a), we plot the macroscopic density at T = 7 . µ ∼ constant = 1, (c) µ ∼ q , (d) µ ∼ q , (e) µ ∼ /q . In (f) we set q constant but maintainanisotropy through k µ = 50 e − . x − +( y − ) and θ µ = π/ µ = q , i.e. situation(62). There is no drift, but anisotropy arises through the directional bias in q . This generatesa conflict in which cells preferentially choose the dominating fibre alignment but, when facingthose directions, turn more frequently. Consequently the anisotropy becomes orthogonalto the dominating fibre alignment. In Fig. 3(d) we consider case (61) by choosing µ = √ q . Again, preferential movement along an axis is countered by increased turning, but theweighting now generates isotropic diffusion. Further, there is a non-trivial drift E T . Thedynamics are similar to those in Figure 3(c), yet the pattern is more symmetric due to theisotropic diffusion. In Fig. 3(e) we consider case (63), i.e. where µ = 1 /q . Cells now turn lessfrequently along the directions of preferential alignment and, intuitively, we should expectenhanced anisotropic diffusion along certain axes. Simulations confirm this, with the cellsremaining even more tightly aligned along the dominating fibre orientation as compared toFig. 3(b). Finally, in Fig. 3(f) the fibre network does not impact on orientation, but cellsoriented along the axis ( π/ , π/
4) turn more frequently. Diffusion remains anisotropic aspredicted by (60).
Test 2 was designed to show the taxis induced by the spatial variability of the turningrate µ and simulation results are reported in Fig. 3 (c). Specifically, we shift the peak of theturning rate away from the centre of the domain, to the point represented by the green dot.There is a subsequent tendency of cells to avoid this new location, coherent with equation25 a) (b) (c)(d) (e) (f) Figure 3:
Test 1 . Anisotropic/isotropic spread according to the choices of q and µ . See textfor details. (a) Initial macroscopic density ¯ p (0 , x ), (b-f) macroscopic density ¯ p ( t, x ) at t = 7 . p ( t, x ) = 0 .
1. In (b-e) q is as described inthe text and the turning rate is given by: (b) µ = 1, (c) µ = q , (d) µ = √ q , (e) µ = 1 /q . In(f) q = 1 / π while µ is given by (64), with θ µ = π/ k µ = 50 e − . x − +( y − ) . (a) t = 2 . t = 5 (c) t = 7 . Figure 4:
Test 2.
Taxis induced through variable turning rates. q is as described in the leftof Fig. 2 (see text for details), while the peak of the turning frequency µ is shifted to centreon the green dot. Simulations plot the macroscopic density at successive times t = 2 . , , . Test 3 extends these analyses to a more complicated environment, see Fig. 2. Fig. 5(a)plots the initial (macroscopic) cell distribution for the experiments (b) and (c), while Fig.5(d) shows the distribution used for (e) and (f). In Fig. 5(b), we set µ ∼ q so that cellsorient and migrate along fibres but this is counterbalanced by turning more frequently when26 a) t = 0 (b) t = 7 . t = 15(d) t = 0 (e) t = 7 . t = 7 . Figure 5:
Test 3 . Dynamics for q given by the cross structure in Fig. 2 (left). (a) Initialmacroscopic density (centered at (5 , , q is given by(65), see text for details, while µ has a similar structure, but for (b,e) θ µ = θ q and for (c,f) θ µ = θ ⊥ q .moving in those directions. Spread is subsequently inhibited by the cross arrangement. InFig. 5(c) we consider instead θ µ = θ ⊥ q , so that particles both follow the fibres and turn lessfrequently when moving in their direction. As expected, there is a clear tendency of cellsto follow the cross structure. The second row performs the same simulations, but relocatingthe initial cell distribution (now centred at (4 , k µ to (3 , µ experienced by cells. Fig. 5 (e) shows even more clearly the inhibition resulting from fibresat the centre of the cross, while in (f) cells are seen to spread rapidly along the arms once ithas been reached. Note the different time scales between Fig. 5(c) and Fig. 5(f). As an application-oriented investigation we return to the cell migration studies of Doyle etal [14], illustrated in Figure 1. These experiments rely on a “photopatterning” technique,allowing fabrication of imprinted fibronectin micro-structures. Laying down parallel alignedstripes imitates aligned fibres and, confronted by such environments, migratory cells (fibrob-lasts and keratinocytes) orient accordingly, extending protrusions and forming the adhesiveattachments that allows movement along the alignment axis. Significantly, highly alignedenvironments lead to a substantial increase in velocity, for example two-fold (for fibroblasts)27r three-fold (for keratinocytes) over corresponding movements on an unaligned surface.While the authors of [14] do not explicitly measure the turning rate µ ( x, w ), they domeasure net velocity and persistence. We note that the emphasis of the experiments in [14]was on cell orientation and morphology, not so much turning rates, so available data remainssparse. Subsequently, rather than a detailed attempt of model fitting, we perform a qualitativecomparison to explore how direction-dependent turning rates will impact on the movementpatterns of cells on different surfaces.We specifically focus on the “transition” experiment, schematically illustrated in Figure1. Here quasi-1D regions were interrupted by two isotropic 2D forms: Case A , a uniformly isotropic region, or
Case B , a region of perpendicular and criss-crossing stripes. As indicatedearlier, cells that enter the isotropic region from a neighbouring quasi-1D region round-up andextend protrusions in multiple directions. For
Case A the cell’s net movement is dramaticallyreduced, losing its direction and subsequently performing what appears as an unbiased randomwalk. In
Case B movement is also arrested and reorientation occurs, but there can be asubsequent significant movement along one of the two perpendicular directions, followed byfurther reorientations. Overall, the translocations of the cell in
Case B seem to be significantlylonger. Here we will show that a direction-dependent turning rate can generate this distinctbehaviour.We adopt a two-pronged approach for the analysis, computing first the macroscopic diffu-sion tensor for cells migrating in the completely unaligned tissue of case A or the criss-crosspattern of case B. While this is a macroscopic-level analysis (and the underlying experimentsare mesoscopic) it will provide valuable evidence of variation in critical movement charac-teristics according to cell orientation/turning behaviour. We then provide simulations of themesoscopic-level transport equation, indicating whether the model can indeed recapitulatethe observations. Note that for convenience we assume an a priori rescaling that fixes thecell speed s , i.e. V = s S . As a control consider a constant turning rate µ and, in turn, a constant mean travel time τ .Here the macroscopic diffusion tensor is computed from formula (58) as D = 1 µ V q . Under
Case A the middle region is completely non-oriented, hence q = π (the uniformdistribution). Then, V q = (cid:90) V w ⊗ w π dw = 2 π I π = 12 I , (66)where I denotes the identity matrix. The criss-cross configuration of Case B can be describedby combining two bi-modal von-Mises distributions, in perpendicular directions e = (1 , e = (0 ,
1) but with the same concentration parameter k : q ( w ) = 12 ( bvM k,e + bvM k,e ) . (67)Following standard calculations (e.g. see [23]), we obtain V q = 12 (cid:18) − I ( k ) I ( k ) (cid:19) I + I ( k ) I ( k ) 12 (cid:0) e e T + e e T (cid:1) . e e T = (cid:18) (cid:19) , and e e T = (cid:18) (cid:19) , so V q = 12 I − I ( k )2 I ( k ) I + I ( k )2 I ( k ) I = 12 I , (68)which coincides exactly with the calculation for Case A , i.e. (66). Therefore, under a constantturning rate there should be no macroscopic difference between Case A and Case B, even ifcells bias their orientation along the criss-crossed fibres.
We now consider w -dependence in the turning rate, i.e. µ ( w ). Specifically, we choose µ suchthat the rate of turning is reduced if cells are migrating along the direction of dominatingalignment. For the analysis we focus only on the central region, so both Case A and
Case B can be regarded as spatially homogeneous (i.e. not depending on x ) and we can use equation(59): D = (cid:90) w ⊗ w Tµ dw, T = qCµ , C ( x ) = (cid:90) V qµ dw. We again choose q to be a combination of bimodal von-Mises distributions in the two per-pendicular directions e and e , as given in (67). However, now we assume that µ ∼ q − , sothat µ = 12 π bvM k,e + bvM k,e ) . where we chose the normalization constant to be (2 π ) − such that in the isotropic limit of k → k → µ = 12 π π + π = 1 . Then Tµ = π C ( bvM k,e + bvM k,e ) . (69)For this choice of µ and q we can compute the normalization constant C as C ( x ) = (cid:90) V qµ dw = (cid:90) V π bvM k,e + bvM k,e ) dw = π πI ( k )) (cid:90) V πI (2 k )( bvM k,e + bvM k,e ) + 4+8 πI ( √ k ) (cid:18) bvM √ k, e e √ + bvM √ k, e − e √ (cid:19) dw = 14 I ( k ) (cid:0) I (2 k ) + 1 + 2 I ( √ k ) (cid:1) . (70)We note the isotropic limit lim k → C ( x ) = 1, since I (0) = 1.29ext we compute the third power in (69) Tµ = π C (cid:18) I (3 k ) I ( k ) (cid:104) bvM k,e + bvM k,e (cid:105) + 3 (cid:104) bvM k,e + bvM k,e (cid:105) +3 (cid:104) bvM k, e + e + bvM k, e − e + bvM k,e +2 e + bvM k,e − e +2( bvM k,e + bvM k,e ) (cid:105)(cid:17) Notably, vectors 2 e + e etc are not unit vectors so we rescale: ζ = 1 √ , T , ζ = 1 √ , − T , ξ = 1 √ , − T , ξ = 1 √ , T . For the example 2 e + e this gives bvM k, e + e = I ( √ k ) I ( k ) bvM √ k,ζ , and similar for the other terms. With unit vectors everywhere which, moreover, are pairwiseperpendicular ( ζ · ζ = 0 , ξ · ξ = 0) we have Tµ = π C (cid:18) I (3 k ) I ( k ) (cid:104) bvM k,e + bvM k,e (cid:105) + 3 (cid:104) bvM k,e + bvM k,e (cid:105) +9( bvM k,e + bvM k,e ) + 3 I ( √ k ) I ( k ) (cid:16) bvM √ k,ζ + bvM √ k,ζ (cid:17) +3 I ( √ k ) I ( k ) (cid:16) bvM √ k,ξ + bvM √ k,ξ (cid:17)(cid:17) . Noting that the second moment of bimodal von Mises distributions with pairwise perpendic-ular unit vectors is I (see (68)), we find a diffusion tensor D = d ( k ) I , d ( k ) = π C (4 πI ( k )) (cid:32) I (3 k ) I ( k ) + 6 I ( √ k ) I ( k ) + 9 (cid:33) . Substituting the normalisation constant C from (70) we obtain d ( k ) = I (3 k ) + 9 I ( k ) + 6 I ( √ k )8 I ( k )( I (2 k ) + 2 I ( √ k ) + 1) . Again we consider the isotropic limitlim k → d ( k ) = 1 + 9 + 68(1 + 2 + 1) = 12 , which has the correct scaling as for the isotropic case. The diffusion coefficient d ( k ) is plottedin Figure 6, where we observe that for small k there is negligible change but for k > d ( k ) as function of k for k ∈ [0 ,
10] on the left and k ∈ [0 , k , but then grows morerapidly for k (cid:38)
3. For large k the increase is approximately ∝ √ k . The above analysis indicates that a direction-dependent turning rate coupled to a criss-crossfibre network substantially increases the diffusion, compared to either a constant turningrate or a completely isotropic network. To simulate this situation, we consider a rectangulardomain Ω = [0 , × [0 ,
5] with V = S ( s = 1) and define local environments correspondingto those illustrated in Figure 1. Specifically, to describe the parallel alignment in the leftand right regions, we set q ( x, w ) = e kw · θ / πI ( k ) with k = 25 and µ ( x, w ) = 1 when x < x >
10. To replicate the completely isotropic scenario of
Case A we define the centralregion (5 ≤ x ≤
10) by q ( x, w ) = 1 / π while for the criss-cross network of Case B we choose q ( x, w ) = ( bvM k,e + bvM k,e ) / k tobe a variable model parameter. Consequently, for Case A µ ( w ) = 1 for the full domain, whilein Case B we chose µ ( x, w ) = 1 / (2 πq ) if 5 ≤ x ≤
10. We initialise the density distribution asuniformly distributed in S with the initial macroscopic density (¯ p ) as defined in 3(a), butcentered at the coordinate (3 . , . k = 0 , , ,
25, respectively.Under Case A (Figure 7 first row), as cells reach the isotropic region we observe a gradualdiffusive-like spread, consistent with the earlier analysis. Under Case B, but setting k =0 (Figure 7 second row), generates equivalent behaviour: in line with the prediction thatisotropic criss-cross networks do not alter the macroscopic dynamics when the turning rate isconstant . For k >
0, however, we observe distinct behaviour. This is minimal for small k (e.g. k = 3, Figure 7 third row) but becomes clear for large k , e.g. k = 10 (fourth row) and k = 25(fifth row), respectively. A noteworthy phenomenon lies in the “droplet” detaching from themain swarm for high anisotropy parameter values ( k = 10 , q ( x, w ) = e kw · θ / (2 πI ( k )) with k = 25 and µ ( x, w ) = 1 for x < x >
10. Firstrow:
Case A with q ( x, w ) = 1 / (2 π ) if 5 ≤ x ≤
10 and µ ( w ) = 1 on the full domain. Secondto final row: q ( x, w ) = ( bvM k,e + bvM k,e ) / µ ( x, w ) = 1 / (2 πq ) if 5 ≤ x ≤
10, so that s = 1. In particular: (second row) k = 0, (third row) k = 3, (fourth row) k = 10, (final row) k = 25. 32 Conclusions
In this paper we have analysed a transport equation for cell migration along oriented fibres,extending the model proposed in [22] to include a turning rate that depends on the microscopicvelocity of the cells and, in turn, on the anisotropic structure of their environment. This keyextension admits a more nuanced and realistic description for how a migratory cell populationresponds to alignment of the environment, with several recent studies investigating the impactof oriented collagen fibre networks on the orientation, speed and persistence of movement ofcells (e.g. see [42, 41, 40]).Formally, the resulting equation is a non-homogeneous linear Boltzmann equation witha micro-reversible process in which the cross-section is factorised into the distribution of thefibres and the direction-dependent turning rate. The dependence of the turning rate on theorientation alters the entire mathematical set up with respect to [22]. First, the equilibriumdistributions of the system now depend on both the turning rate and on the transition prob-ability. Consequently, the null-space of the turning operator needs to be equipped with adirection-dependent integration weight. Moreover, the average of the equilibrium state, de-fined in Eq. (11), becomes linked to both the fibre distribution and the turning frequency.This implies new solvability conditions for the parabolic and hyperbolic limits. We studymany special cases, some of which are relevant for applications while others provide a theo-retical connection to previous results.Further, the orientation-dependent turning rate permits the definition of a new adjointpersistence (13), taking into account the cell persistence encoded in the turning rate. Con-sequently, direction-dependent turning rates are capable of generating a persistent randomwalk even under a symmetric configuration of fibres. In this framework, there is no direc-tional persistence if (11) vanishes. This holds true, for example, if Eq. (19) holds true, whichmeans that fibres are bi-directed and cells having a certain direction have the same turningfrequency regardless of the sense in which they are travelling on that direction.To illustrate the broader relevance of the extended framework, we considered an appli-cation to cell migration across precisely engineered network arrangements, for example asfabricated in the studies of [14]. Under the the original framework of [22] the two configura-tions shown in Figure 1 generate equivalent behaviour, yet extending to direction-dependentturning rates could result in a markedly different response. Specifically, under the criss-crossnetwork a markedly faster passage could be observed which, translated to the macroscopiclevel, yielded enhanced diffusion. Clearly, such behaviour has potential to significantly alterthe predictions from modelling invasion pathways in complex anisotropic tissues, for exampleglioma invasion in the central nervous tissue [35, 16, 48].Within the current work we have restricted to movement under negligible modificationof the network, as in amoeboid movement or cell migration on engineered fibronectin strips.Under in vivo mesenchymal migration, however, contact-guided movement can be coupledto significant matrix remodelling, for example fibres becoming aligned along the migratorypath; simulations of the simpler transport model in this scenario revealed symmetry breakingbehaviour, with cells forming and migrating along a network of aligned “cellular highways”[22, 34]. The addition of direction-dependent turning has clear potential to impact on suchpattern formation scenarios and a key aim of future studies will be to investigate this phe-nomenon.
Acknowledgements:
This work was stimulated by the PhD thesis of Amanda Swan [47], who tarted looking into the parabolic scaling of a velocity dependent turning rate. TH is grateful tosupport through the Natural Sciences and Engineering Research Council of Canada. NL is recipientof a Post-Doc grant of the Italian National Institute of High Mathematics (INdAM) and ackowledgessupport by the Italian Ministry for Education, University and Research (MIUR) through the “Diparti-menti di Eccellenza” Programme (2018- 2022), Department of Mathematical Sciences, G. L. Lagrange,Politecnico di Torino (CUP: E11G18000350001), and support from “Compagnia di San Paolo” (Torino,Italy). References [1] J. Adler. Chemotaxis in bacteria.
Science , 153(3737):708–116, 1966.[2] B. Alberts, A. D. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter.
Molecular Biology of the Cell . Garland Science, Taylor and Francis Group,, 2014.[3] M. Bisi, J. A. Carrillo, and B. Lods. Equilibrium solution to the inelastic Boltzmannequation driven by a particle bath.
Journal of Statistical Physics , 133(5):841–870, 2008.[4] E. O. Budrene and H. C. Berg. Complex patterns formed by motile cells of escherichiacoli.
Nature , 349(6310):630–633, 1991.[5] C. Cercignani.
The Boltzmann Equation and its Applications . Springer, New York, 1987.[6] F. A. C. C. Chalub, P. A. Markowich, B. Perthame, and C. Schmeiser. Kinetic models forchemotaxis and their drift-diffusion limits.
Monatshefte f¨ur Mathematik , 142(1):123–141,Jun 2004.[7] S. Chapman. On the brownian displacements and thermal diffusion of grains suspendedin a non-uniform fluid.
Proceedings of the Royal Society London A , 119:34–54, 1928.[8] E. Cho and Y. J. Kim. Starvation driven diffusion as a survival strategy of biologicalorganisms.
Bulletin of Mathematical Biology , 75:845–870, 2013.[9] J. Chung, Y. J. Kim, O. Kwong, and C. W. Yoon. Biological advection and cross diffusionwith parameter regimes.
AIMS Mathematics , 4(6), 2020.[10] J. C. Dallon, J. A. Sherratt, and P. K. Maini. Mathematical modelling of extracellularmatrix dynamics using discrete cells: fiber orientation and tissue regeneration.
Journalof Theoretical Biology , 199(4):449–471, 1999.[11] P. Degond, T. Goudon, and F. Poupaud. Diffusion limit for non homogeneous and non-micro-reversible processes.
Indiana University Mathematics Journal , 49(3):1175–1198,2000.[12] R. B. Dickinson. A generalized transport model for biased cell migration in an anisotropicenvironment.
Journal of Mathematical Biology , 40(2):97–135, Feb 2000.[13] R. B. Dickinson and R. T. Tranquillo. Stochastic model of biased cell migration based onbinding fluctuations of adhesion receptors.
Journal of Mathematical Biology , 19:563–600,1991. 3414] A. D. Doyle, F. W. Wang, K. Matsumoto, and K. M. Yamada. One-dimensional to-pography underlies three-dimensional fibrillar cell migration.
Journal of Cell Biology ,184(4):481–490, 2009.[15] G.A. Dunn and J.P. Heath. A new hypothesis of contact guidance in tissue cells.
Exper-imental Cell Research , 101(1):1–14, 1976.[16] C. Engwer, Thomas Hillen, M. Knappitsch, and C. Surulescu. Glioma follow white mattertracts: a multiscale DTI-based model.
Journal of Mathematical Biology , 71(3):551–582,2015.[17] F. Filbet and C Yang. Numerical simulation of a kinetic model for chemotaxis.
Kineticand Related Models , 3:B348–B366, 2010.[18] A. Giese, L. Kluwe, B. Laube, H. Meissner, M. E. Berens, and M. Westphal. Migrationof human glioma cells on myelin.
Neurosurgery , 38(4):755–764, 1996.[19] A. Giese and M. Westphal. Glioma invasion in the central nervous system.
Neurosurgery ,39(2):235–252, 1996.[20] T. Goudon and A. Mellet. Homogenization and diffusion asymptotics of the linear Boltz-mann equation.
ESAIM Control Optimisation and Calculus of Variations , 9:371–398, 042003.[21] I. Hecht, Y. Bar-El, F. Balmer, S. Natan, I. Tsarfaty, F. Schweitzer, and E. Ben-Jacob.Tumor invasion optimization by mesenchymal-amoeboid heterogeneity.
Scientific Re-ports , 5(10622), 2015.[22] T. Hillen. M5 mesoscopic and macroscopic models for mesenchymal motion.
Journal ofMathematical Biology , 53(4):585–616, 2006.[23] T. Hillen, A. Murtha, K.J. Painter, and A. Swan. Moments of the von Mises and Fisherdistributions and applications.
Mathematical Biosciences and Engineering , 14(3):673–694, 2017.[24] T. Hillen and H. G. Othmer. The diffusion limit of transport equations derived fromvelocity-jump processes.
SIAM Journal of Applied Mathematics , 61:751–775, 2000.[25] T. Hillen and K. J. Painter. A user’s guide to PDE models for chemotaxis.
Journal ofMathematical Biology , 58(1):183–217, 2008.[26] T. Hillen and K. J. Painter.
Transport and anisotropic diffusion models for movementin oriented habitats , volume 2071, pages 177–222. Springer-Verlag, 2013.[27] B. Lods. Semigroup generation properties of streaming operators with noncontractiveboundary conditions.
Mathematical and Computer Modelling , 42:1441–1462, 2005.[28] N. Loy and L. Preziosi. Kinetic models with non-local sensing determining cell polariza-tion and speed according to independent cues.
Journal of Mathematical Biology , pages1–49, 2019.[29] F. Lutscher and T. Hillen. Homogenization of correlated random walks in heterogeneouslandscapes.
AIMS Mathemtics , 2021. to appear.3530] S. McDougall, J. Dallon, J. Sherratt, and P. Maini. Fibroblast migration and collagendeposition during dermal wound healing: mathematical modelling and clinical implica-tions.
Philosophical Transactions of the Royal Society A: Mathematical, Physical andEngineering Sciences , 364(1843):1385–1405, 2006.[31] H. Othmer and A. Stevens. Aggregation, blowup, and collapse: The ABC’s of taxis inreinforced random walks.
SIAM Journal on Applied Mathematics , 57:1044–1081, 2001.[32] H. G. Othmer, S. R. Dunbar, and W. Alt. Models of dispersal in biological systems.
Journal of Mathematical Biology , 26(3):263–298, 1988.[33] H. G. Othmer and T. Hillen. The diffusion limit of transport equations II: Chemotaxisequations.
SIAM Journal of Applied Mathematics , 62:1222–1250, 2002.[34] K. J. Painter. Modelling cell migration strategies in the extracellular matrix.
Journal ofMathematical Biology , 58(4):511–543, 2008.[35] K. J. Painter and T. Hillen. Mathematical modelling of glioma growth: the use of diffu-sion tensor imaging (DTI) data to predict the anisotropic pathways of cancer invasion.
Journal of Theoretical Biology , 323:25–39, 2013.[36] Kevin J. Painter and Thomas Hillen.
From Random Walks to Fully Anisotropic Dif-fusion Models for Cell and Animal Movement , pages 103–141. Springer InternationalPublishing, Cham, 2018.[37] R. Petterson. Existence Theorems for the Linear, Space-inhomogeneous Transport Equa-tion.
IMA Journal of Applied Mathematics , 30(1):81–105, 1983.[38] R. Pettersson. On solutions to the linear boltzmann equation for granular gases.
Trans-port Theory and Statistical Physics , 33(5-7):527–543, 2004.[39] R. G. Plaza. Derivation of a bacterial nutrient-taxis system with doubly degeneratecross-diffusion as the parabolic limit of a velocity-jump process.
Journal of MathematicalBiology , 78:1681–1711, 2019.[40] A. Ray, R. K. Morford, N. Ghaderi, D. J. Odde, and P. P. Provenzano. Dynamics of 3dcarcinoma cell invasion into aligned collagen.
Integrative Biology , 10(2):100–112, 2018.[41] A. Ray, Z. M. Slama, R. K. Morford, S. A. Madden, and P. P. Provenzano. Enhanceddirectional migration of cancer stem cells in 3d aligned collagen matrices.
BiophysicalJournal , 112(5):1023–1036, 2017.[42] K. M. Riching, B. L. Cox, M. R. Salick, C. Pehlke, A. S. Riching, S. M. Ponik, B. R.Bass, W. C. Crone, Y. Jiang, A. M. Weaver, K.W. Eliceiri, and P. J. Keely. 3d colla-gen alignment limits protrusions to enhance breast cancer cell persistence.
BiophysicalJournal , 107(11):2546–2558, 2014.[43] D. K. Schl¨uter, I. Ramis-Conde, and M. A. J. Chaplain. Computational modeling ofsingle-cell migration: the leading role of extracellular matrix fibers.
Biophysical Journal ,103(6):1141–1151, 2012. 3644] M. Scianna and L. Preziosi. Modeling the influence of nucleus elasticity on cell invasionin fiber networks and microchannels.
Journal of Theoretical Biology , 317:394–406, 2013.[45] M. Scianna and L. Preziosi. A cellular Potts model for the MMP-dependent and-independent cancer cell migration in matrix microtracks of different dimensions.
Com-putational Mechanics , 53:485–497, 2014.[46] M. Scianna, L. Preziosi, and K. Wolf. A Cellular Potts Model simulating cell migrationon and in matrix environments.
Mathematical Biosciences and Engineering , 10:235–261,2013.[47] A. Swan.
An Anisotropic Diffusion Model for Brain Tumour Spread . PhD thesis, Uni-versity of Alberta, 2016.[48] A. Swan, T. Hillen, J. Bowman, and A. Murtha. An anisotropic model for glioma spread.
Bulletin of Mathematical Biology , 80(5):1259–1291, 2017.[49] K. Talkenberger, E. A. Cavalcanti-Adam, A. Voss-B¨ohme, and A. Deutsch. Amoeboid-mesenchymal migration plasticity promotes invasion only in complex heterogeneous mi-croenvironments.
Scientific Reports , 7(9237), 2017.[50] V. Te Boekhorst, L. Preziosi, and P. Friedl. Plasticity of cell migration in vivo and insilico.
Annual Review of Cell and Developmental Biology , 32:491–526, 2016.[51] M. Th´ery. Micropatterning as a tool to decipher cell morphogenesis and functions.
Journal of Cell Science , 123(24):4201–4213, 2010.[52] M. Th. Wereide. La diffusion d’une solution dont la concentration et la temperature sontvariables.
Annales de Physique , 2:67–83, 1914.[53] K. Wolf, I. Mazo, I. Leung, K. Engelke, U. von Andria, E.I. Deryngina, A.Y. Stron gin,E.B. Brocker, and P. Friedl. Compensation mechanism in tumor cell migration:mesenchymal-amoeboid transition after blocking of pericellular proteolysis.