Macroscopic limit of a kinetic model describing the switch in T cell migration modes via binary interactions
MMacroscopic limit of a kinetic model describing the switch in T cellmigration modes via binary interactions
Gissell Estrada-Rodriguez ∗ Tommaso Lorenzi † Abstract
Experimental results on the immune response to cancer indicate that activation of cytotoxic T lym-phocytes (CTLs) through interactions with dendritic cells (DCs) can trigger a change in CTL migrationpatterns. In particular, while CTLs in the pre-activation state move in a non-local search pattern, thesearch pattern of activated CTLs is more localised. In this paper, we develop a kinetic model for such aswitch in CTL migration modes. The model is formulated as a coupled system of balance equations forthe one-particle distribution functions of CTLs in the pre-activation state, activated CTLs and DCs. CTLactivation is modelled via binary interactions between CTLs in the pre-activation state and DCs. Moreover,cell motion is represented as a velocity-jump process, with the running time of CTLs in the pre-activationstate following a long-tailed distribution, which is consistent with a L´evy walk, and the running time ofactivated CTLs following a Poisson distribution, which corresponds to Brownian motion. We formally showthat the macroscopic limit of the model comprises a coupled system of balance equations for the cell densitieswhereby activated CTL movement is described via a classical diffusion term, whilst a fractional diffusionterm describes the movement of CTLs in the pre-activation state. The modelling approach presented hereand its possible generalisations are expected to find applications in the study of the immune response tocancer and in other biological contexts in which switch from non-local to localised migration patterns occurs.
The interaction between dendritic cells (DCs) and cytotoxic T lymphocytes (CTLs) plays a pivotal role inthe immune response to cancer. DCs recognise the antigens expressed by cancer cells and present them toCTLs, which then become selectively activated against those antigens [21, 22]. Growing experimental evidenceindicates that activation of CTLs via antigen presentation by DCs can bring about a switch in CTL migrationmodes [3, 13]. In fact, while CTLs in the pre-activation state move in a non-local search pattern, which enablesthem to rapidly scan DCs for the presence of possible tumour antigens, the search pattern of activated CTLsis more localised. This allows activated CTLs to stay within a confined area for longer, thus facilitating theirencounter with tumour cells expressing the antigens they have been activated against.Stochastic individual-based models of immune response to cancer taking explicitly into account this differencein movement between CTLs have recently been developed [14, 15]. In these models, cell motion is describedas a space-jump process [17]. In particular, CTLs in the pre-activation state undergo a space-jump processconsistent with a L´evy walk, whereas a space-jump process corresponding to Brownian motion is used todescribe the movement of activated CTLs. Such individual-based models enable representation of biological ∗ Sorbonne Universit´e, CNRS, Universit´e de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris, France([email protected]) † Department of Mathematical Sciences “G. L. Lagrange”, Dipartimento di Eccellenza 2018-2022, Politecnico di Torino, 10129Torino, Italy ([email protected]) a r X i v : . [ q - b i o . CB ] J u l rocesses at the level of single cells and account for possible stochastic variability in cell dynamics, which allowfor greater adaptability and higher accuracy in mathematical modelling. However, as the numerical explorationof these models requires large computational times for clinically relevant cell numbers and the models are notanalytically tractable, it is desirable to derive corresponding deterministic continuum models in a suitable limit.In this paper, integrating the ideas proposed in [14, 15] with the modelling approach presented in [8, 9], wedevelop a kinetic model for the switch in CTL migration modes that is caused by activation through interactionswith DCs. Cells are grouped into three populations: CTLs in the pre-activation state ( i.e. inactive CTLs),activated CTLs and DCs. In the model, DCs are assumed to present a given tumour antigen on their surface sothat they can activate inactive CTLs by contact. Since the focus of this study is on the mathematical modellingof the change in CTL migration mode upon activation, we do not take into account biological processes involvingcell division and death. Furthermore, for simplicity, we do not consider loss of effector functions leading activatedCTLs to re-enter a pre-activation state [23].The model is formulated as a coupled system of balance equations for the one-particle distribution functions ofthe three cell populations. CTL activation is modelled as a process of population switching among CTLs inducedby binary interactions between inactive CTLs and DCs. Moreover, cell motion is represented as a velocity-jumpprocess [17], with the running time of inactive CTLs following a long-tailed distribution, which is consistent witha L´evy walk [8, 9], and the running time of activated CTLs following a Poisson distribution, which correspondsto Brownian motion. Using a method similar to that previously employed in [8], we formally show that themacroscopic limit of this model comprises a coupled system of balance equations for the cell densities, wherebyactivated CTL movement is described via a classical diffusion term, whilst a fractional diffusion term describesthe movement of CTLs in the pre-activation state.The paper is organised as follows. In Section 2, we introduce the modelling strategies and the main assump-tions that are used to describe the spatio-temporal dynamics of CTLs and DCs at the scale of single cells,which provide a microscopic representation of the biological system. In Section 3, we present the kinetic model,which constitutes a mesoscopic analogue of the underlying microscopic scale model. In Section 4, we derive themacroscopic limit of a suitably rescaled version of the kinetic model. Section 5 concludes the paper providinga brief overview of possible research perspectives. Biological system and cell populations
We label the three cell populations by a letter h ∈ { A, D, I } , thatis, activated CTLs are labelled by h = A , DCs are labelled by h = D and inactive CTLs are labelled by h = I .We denote the total number of cells in the system by N = N D + N T , where N D ∈ N is the number of DCs and N T ∈ N is the total number of CTLs. Moreover, we describe the number of inactive and activated CTLs in thesystem at time t ∈ R + by means of the functions N I ( t ) and N A ( t ), respectively, with N I ( t ) + N A ( t ) = N T forall t . Mathematical representation of individual cells
Every individual cell is modelled as a sphere of diameter (cid:37) ∈ R ∗ + and is labelled by an index i = 1 , . . . , N . The phase-space state of the i th cell is represented by a pair( x i , v i ), where the vector x i ∈ R n describes the position of the centre of the cell and the vector v i ∈ V ⊂ R n ,with V := { v i ∈ R n : | v i | = 1 } ( i.e. V is the unit n -sphere), represents the direction of the cell velocity.Moreover, the magnitude of the cell velocity is assumed to be constant and is denoted by c ∈ R ∗ + . The value of n = 1 , , .1 Description of cell motion Velocity-jump process
We describe the motion of a cell labelled by an index i as a run-and-tumble processwith run time τ i ∈ R ∗ + and running probability ψ i ( x i , τ i ), where 0 < ψ i ( · , · ) ≤ ∂ τ i ψ i ( · , · ) ≤
0. The runningprobability ψ i ( x i , τ i ) correlates with the stopping rate β i ( x i , τ i ) through the relations given by the followingdefinition [8] ψ i ( x i , τ i ) := exp (cid:18)(cid:90) τ i β i ( x i , s ) d s (cid:19) , β i = ϕ i ψ i with ϕ i := − ∂ τ i ψ i . (2.1)Hence, starting at position x i at time t , the i th cell will continue moving along a straight path in the directiongiven by the vector v i with constant speed c for a period of time τ i , after which it may stop at rate β i ( x i , τ i ).The cell will then instantaneously resume moving in a new randomly selected direction given by a vector ¯ v i ,which is prescribed by a turning kernel (cid:96) ( x i , t, v i ; ¯ v i ) – i.e. cells undergo a velocity-jump process [17]. Running probability
The running probability ψ i ( x i , τ i ) determines the distribution of the running time τ i and depends on the way in which the i th cell moves. On the basis of experimental evidence reported in [3, 7],we assume that: inactive CTLs move in a non-local search pattern corresponding to trajectories that arecharacterised by a strong presence of long runs, which enable them to cover larger areas; activated CTLs andDCs move in a more localised search pattern. In particular, building upon the modelling approach presentedin [14, 15], we describe the motion of activated CTLs and DCs as a Brownian motion, whereas we let inactiveCTLs undergo superdiffusive motion consistent with a L´evy walk, whereby the mean square displacement growsnonlinearly with time. In particular, the mean square displacement at time t is proportional to t / α , where α ∈ (1 ,
2) is the L´evy exponent. We recall that α = 1 and α = 2 correspond to ballistic motion and classicaldiffusion, respectively.Under these assumptions, if the i th cell belongs to population A or population D , we let the value of therunning time τ i follow a Poisson distribution [18]. Hence, under the additional simplifying assumption that cellsin populations A and D are characterised by the same stopping rate, which is assumed to be constant and thusmodelled by a parameter b ∈ R ∗ + , we use the following definition of the running probability ψ i ( x i , τ i ) ≡ ψ i ( τ i ) := exp ( − b τ i ) , ϕ i ( x i , τ i ) ≡ ϕ i ( τ i ) := b exp ( − b τ i ) . (2.2)On the other hand, if the i th cell belongs to population I , we let the value of the running time τ i follow along-tailed distribution, and we define the running probability along the lines of [8] as ψ i ( x i , τ i ) := (cid:16) τ ( x i ) τ ( x i ) + τ i (cid:17) α , ϕ i ( x i , τ i ) := α τ ( x i ) α ( τ ( x i ) + τ i ) α +1 , α ∈ (1 , . (2.3)Here, the function τ ( x i ) ≥ Turning kernel and turning operator
We consider the case where the new direction of cell motion givenby ¯ v i is symmetrically distributed with respect to the original direction given by v i and, therefore, we let theturning kernel (cid:96) ( x i , t, v i ; ¯ v i ) satisfy the following assumptions [1] (cid:96) ( x i , t, v i ; ¯ v i ) ≡ (cid:96) ( x i , t, | ¯ v i − v i | ) , (cid:90) V (cid:96) ( · , · , | v i − e | ) d v i = 1 , (2.4)where e = (1 , , . . . , ∈ R n is a unit vector. We remind the reader that we consider DCs presenting a given tumour antigen on their surface.
3e let the integral operator T i be a turning operator such that for all test functions φ ( v i ) T i [ φ ]( · , · , ¯ v i ) = (cid:90) V (cid:96) ( · , · , v i ; ¯ v i ) φ ( v i ) d v i , (2.5)where (cid:96) is the turning kernel defined via (2.4). Since (cid:90) V (cid:96) ( · , · , · ; ¯ v ) d¯ v = 1, we have (cid:90) V ( − T i )[ φ ]( · , · , ¯ v i ) d¯ v i = 0 , (2.6)where is the identity operator. Focusing on a biological scenario of relatively low cell densities ( cf. the scaling and assumptions introducedin Section 4.1), we consider only the effect of binary interactions between cells, thus neglecting interactionsinvolving more than two cells [4, 17, 20]. We assume that interactions between a test cell in the phase-spacestate ( x i , v i ) and a field cell in the phase-space state ( x j , v j ) occur when the field cell is in the domain ofinteraction of the test cell, which is defined as the setΩ j ( x i ) := { x j ∈ R n : | x i − x j | ≥ (cid:37) } ≡ R n \ B (cid:37) ( x i ) , (2.7)where B (cid:37) ( x i ) denotes the ball of radius (cid:37) centred at x i . We model interactions between cells causing a changein the velocity of the test cell as elastic collisions, whereby the post-collision velocity v (cid:48) i of the test cell is givenby v (cid:48) i = v i − v i · ν ) ν with ν := x i − x j | x i − x j | , (2.8)where ν is the normal vector at the point of interaction ( i.e. ν is the unit normal that points outward fromΩ j ( x i ) and inward to B (cid:37) ( x i )) [4].In particular, we consider the binary interactions between test cells and field cells that are summarised bythe schematics in Figure 1, which correspond to the following assumptions. Assumption 1 ( Interactions between inactive CTLs and DCs).
We model activation of CTLs uponinteraction with DCs by assuming that, when a test cell in population I interacts with a field cell in population D , the test cell switches from population I to population A ( i.e. the interaction is destructive for population I and creative for population A ) with probability ζ ∈ (0 ,
1) and its velocity remains unchanged. If activation doesnot occur, event that happens with probability 1 − ζ , the test cell remains in population I ( i.e. the interactionis conservative) and acquires the post-collision velocity defined via (2.8). Assumption 2 ( Interactions between activated CTLs and DCs).
We assume that when a test cell inpopulation A interacts with a field cell in population D , the test cell acquires the post-collision velocity definedvia (2.8) and the interaction is conservative. Assumption 3 ( Interactions between DCs and CTLs).
We assume that if a test cell in population D interacts with a field cell in population I , the test cell will acquire the post-collision velocity defined via (2.8)and the interaction will be conservative. For simplicity, we do not take into account the effect of interactionsbetween test cells of population D and field cells of population A .4igure 1: Schematics of cell-cell interactions corresponding to Assumptions 1-4. Assumption 4 ( Other cell-cell interactions).
Since here we are primarily interested in the mathematicalmodelling of the switch in T cell migration modes mediated by interactions between inactive CTLs and DCs,for simplicity we neglect interactions between cells of the same population and we do not take into account theeffect of interactions between cells of population I and cells of population A . In this section, we derive the mesoscopic scale model corresponding to the microscopic scale description presentedin Section 2, which comprises a system of transport equations for the one-particle distribution functions ofinactive CTLs, activated CTLs and DCs.
The state of the system at time t is described by the N -particle distribution function f N ( x , . . . , x N , t, v , . . . , v N , τ , . . . , τ N )[4, 20]. In the case where cell dynamics at the microscopic scale obey the rules presented in Section 2, the evo-lution of f N is governed by the following transport equation [12] ∂ t f N + N (cid:88) i =1 (cid:16) ∂ τ i f N + c v i · ∇ x i f N (cid:17) = − N (cid:88) i =1 β i f N (3.1)posed on Ω N × R ∗ + × V N × R ∗ N + , withΩ N := { ( x , ..., x N ) ∈ R n × N : | x i − x j | ≥ (cid:37) ∀ i, j } . We consider the transport equation (3.1) subject to smooth, compactly supported initial conditions at t =0, specular reflective boundary conditions corresponding to elastic collisions on ∂ Ω N , and suitable Dirichletboundary conditions at τ i = 0 linked to the running probability ψ i for i = 1 , . . . , N . In the mathematical5ramework given by (3.1), the probability of finding at position x and at time t the cell labelled by the index1 that is moving in direction v for a period of time τ is related to the one-particle marginal f ( x , t, v , τ ) = 1 | V | N − (cid:90) [0 ,t ] N − (cid:90) Ω N − ( x ) (cid:90) V N − f N ( x , . . . , x N , t, v , . . . , v N , τ , . . . , τ N ) × d v d x d τ . . . d v N d x N d τ N . Here, | V | is the surface area of the unit sphere V and Ω N − ( x ) := { ( x , . . . , x N ) ∈ R n × N − : ( x , x , . . . , x N ) ∈ Ω N } .A comprehensive description of cell dynamics would in principle require considering possible interactionsbetween all cells. However, as mentioned earlier, we focus on a biological scenario whereby cell densities arerelatively low, as per the scaling and assumptions introduced in Section 4.1. In such a low-density regime,interactions between more than two cells at a time can be neglected [4, 17, 20]. Therefore, we truncate thehierarchy of equations corresponding to (3.1) at the second order by integrating out cells 3 , . . . , N from the N -particle distribution function f N ( x , . . . , x N , t, v , . . . , v N , τ , . . . , τ N ). Two-particle distribution functions
We denote by f hk ( x h , x k , t, v h , v k , τ h , τ h ) with h, k ∈ { A, D, I } and k (cid:54) = h the two-particle distribution function associated with:- a test cell of population h in the generic phase-space state ( x h , v h ) ∈ R n × V, with generic run time τ h ∈ [0 , t ] and stopping rate β h ( x h , τ h ) defined via (2.1);- a field cell of population k in the generic phase-space state ( x k , v k ) ∈ R n × V, with generic run time τ k ∈ [0 , t ] and stopping rate β k ( x k , τ k ) defined via (2.1).Truncating the hierarchy of equations corresponding to (3.1) at the second order, we obtain the followingtransport equation for f hk ( x h , x k , t, v h , v k , τ h , τ h )( ∂ t + ∂ τ h + ∂ τ k + c v h · ∇ x h + c v k · ∇ x k ) f hk = − ( β h + β k ) f hk (3.2)posed on Ω × R + × V × R ∗ , withΩ := { ( x h , x k ) ∈ R n × : | x h − x k | ≥ (cid:37) ∀ h, k } . (3.3)This equation is subject to a smooth, compactly supported initial condition at t = 0, specular reflective boundaryconditions corresponding to elastic collisions on ∂ Ω , and Dirichlet boundary conditions at τ h = 0 and τ k = 0linked to the running probabilities ψ h and ψ k , respectively. One-particle distribution functions
Given the two-particle distribution function˜˜ f hk ( x h , x k , t, v h , v k ) := (cid:90) t (cid:90) t f hk d τ h d τ k , (3.4)the one-particle distribution function of population h is given by p h ( x h , t, v h ) := 1 | V | (cid:90) Ω k ( x h ) (cid:90) V ˜˜ f hk d v k d x k , (3.5)with the set Ω k ( x h ) defined via (2.7). Moreover, we will consider the weighted two-particle distribution functiongiven by ˜˜ f β z hk ( x h , x k , t, v h , v k ) := (cid:90) t (cid:90) t β z f hk d τ h d τ k , z ∈ { h, k } , (3.6)6nd the weighted one-particle distribution function given by p β h h ( x h , t, v h ) := 1 | V | (cid:90) Ω k ( x h ) (cid:90) V ˜˜ f β h hk d v k d x k . (3.7) Transport equations for two-particle distribution functions
The dynamics of the two-particle distri-bution functions f ID , f AD and f DI are governed by the following specific forms of transport equation (3.2),( ∂ t + ∂ τ I + ∂ τ D + c v I · ∇ x I + c v D · ∇ x D ) f ID = − ( β I + β D ) f ID , (3.8)( ∂ t + ∂ τ A + ∂ τ D + c v A · ∇ x A + c v D · ∇ x D ) f AD = − ( β A + β D ) f AD , (3.9)( ∂ t + ∂ τ D + ∂ τ I + c v D · ∇ x D + c v I · ∇ x I ) f DI = − ( β I + β D ) f DI , (3.10)which are posed on Ω × R ∗ + × V × R ∗ . Starting from transport equations (3.8)-(3.10) and using the methodemployed in [8], it is possible to show (see Appendix A) that the two-particle distribution functions ˜˜ f ID , ˜˜ f AD and ˜˜ f DI given by (3.4) satisfy the following transport equations( ∂ t + c v I · ∇ x I + c v D · ∇ x D ) ˜˜ f ID = − ( − T I )[ ˜˜ f β I ID ] − ( − T D )[ ˜˜ f β D ID ] , (3.11)( ∂ t + c v A · ∇ x A + c v D · ∇ x D ) ˜˜ f AD = − ( − T A )[ ˜˜ f β A AD ] − ( − T D )[ ˜˜ f β D AD ] (3.12)and ( ∂ t + c v D · ∇ x D + c v I · ∇ x I ) ˜˜ f DI = − ( − T D )[ ˜˜ f β D DI ] − ( − T I )[ ˜˜ f β I DI ] , (3.13)posed on Ω × R ∗ + × V . Here, T I , T D and T A are the turning operators defined via (2.5), and ˜˜ f β h hk and ˜˜ f β k hk arethe weighted two-particle distribution functions given by (3.6). Transport equation for p h Starting from transport equation (3.2) and building upon the method presentedin [8], it is possible to show (see Appendix B) that the one-particle distribution function p h ( x h , t, v h ) givenby (3.5) satisfies the following transport equation ∂ t p h + c v h · ∇ x h p h = − ( − T h )[ p β h h ] + Q hk , x h ∈ R n , t ∈ R + , v h ∈ V . (3.14)Here, the turning operator T h is defined via (2.5), the weighted one-particle distribution function p β h h ( x h , t, v h )is given by (3.7) and Q hk ( x h , t, v h ) := c | V | (cid:90) ∂ B (cid:37) ( x h ) (cid:90) V ν · ( v h − v k ) ˜˜ f hk d v k d σ . (3.15)In (3.15), ν is the unit normal defined in (2.8) and d σ denotes the surface element.The first term on the right-hand side of transport equation (3.14) represents the rate of change of the one-particle distribution function due to cell movement, while the term Q hk is the rate of change due to interactionsbetween cells. The specific forms of these terms depend, respectively, on the way in which cells move and theinteractions they undergo, as discussed in the remainder of this section.7 xpressions for p β h h The specific form of the first term on the right-hand side of transport equation (3.14)depends on the expression for p β h h which, in turn, will depend on the definition of the stopping rate β h .When cells move in a local search pattern ( i.e. for h = A and h = D ), the stopping rate β h is defined via (2.1)and (2.2). In this case, inserting the definition of β h into (3.7) yields p β h h ( x h , t, v h ) = b p h ( x h , t, v h ) . (3.16)On the other hand, when cells move in a non-local search pattern ( i.e. for h = I ), the stopping rate β h isdefined via (2.1) and (2.3). In this case, it is possible to show (see Appendix C) that p β h h ( x h , t, v h ) = B [ p h ]( x h , t, v h ) , (3.17)where B is a convolution operator such that B [ p h ]( x h , t, v h ) = (cid:90) t B ( x h , t − s ) p ( x h − ( c v h + b )( t − s ) , s, v h ) d s , (3.18)with B being defined through its Laplace transform in time ˆ B asˆ B ( x h , λ + b + c v h · ∇ x h ) = ˆ ϕ h ( x h , λ + b + c v h · ∇ x h )ˆ ψ h ( x h , λ + b + c v h · ∇ x h ) . (3.19)Here, λ is the Laplace variable, ˆ ϕ h and ˆ ψ h are the Laplace transforms in τ h of the functions ϕ h and ψ h definedvia (2.3), and the parameter b is defined via (2.2). Expressions for Q hk Following [8, 10], we first note that when a test cell in the phase-space state ( x h , v h )interacts with a field cell in the phase-space state ( x k , v k ) we have | x h − x k | = (cid:37) . Hence, the normal vector at thepoint of physical contact between the interacting cells, ν ∈ V, defined via (2.8) can be written as ν = ( x h − x k ) /(cid:37) ,that is, x k = x h − ν(cid:37) . As a result, using the fact that B (cid:37) = (cid:37) V along with the change of variable ν (cid:55)→ − ν , werewrite (3.15) as Q hk ( x h , t, v h ) := c | V | (cid:90) ∂ B (cid:37) ( x h ) (cid:90) V ν · ( v h − v k ) ˜˜ f hk d v k d σ = − c | V | (cid:37) n − (cid:90) V (cid:90) V ν · ( v h − v k ) ˜˜ f hk ( x h , x h + ν(cid:37), t, v h , v k ) d v k d ν . (3.20)Following [8], we also note that V ≡ V + hk ∪ V − hk withV + hk := { ν ∈ V : ν · ( v h − v k ) > } , V − hk := { ν ∈ V : ν · ( v h − v k ) < } ≡ {− ν ∈ V : ν · ( v h − v k ) > } . (3.21)Therefore, a test cell moving in direction v h and a field cell moving in direction v k will move toward each otherif ν ∈ V + and away from each other if ν ∈ V − .Under Assumptions 1-4, denoting the post-collision directions corresponding to v h and v k by v (cid:48) h and v (cid:48) k ,which are defined via (2.8), we consider two different types of interactions between cells:- conservatives interactions, between a test cell of population h and a field cell of population k , wherebyboth cells remain in their original populations upon interaction and acquire the post-collision velocities;8 population-switching interactions, between a test cell in population h and a field cell in population k ,whereby the test cell switches from its original population to a different one upon interaction and the cellvelocities remain unchanged.From (3.20), we define the rate of change of the one-particle distribution function p h ( x h , t, v h ) due to conser-vative interactions as K hk ( x h , t, v h ) := − c | V | (cid:37) n − N k ( t ) (cid:104) (cid:90) V + hk (cid:90) V ν · ( v h − v k ) ˜˜ f hk ( x h , x h + ν(cid:37), t, v h , v k ) d v k d ν + (cid:90) V − hk (cid:90) V ν · ( v h − v k ) ˜˜ f hk ( x h , x h + ν(cid:37), t, v h , v k ) d v k d ν (cid:105) = c | V | (cid:37) n − N k ( t ) (cid:90) V + hk (cid:90) V ν · ( v h − v k ) (cid:104) ˜˜ f hk ( x h , x h − ν(cid:37), t, v (cid:48) h , v (cid:48) k ) − ˜˜ f hk ( x h , x h + ν(cid:37), t, v h , v k ) (cid:105) d v k d ν , (3.22)with N k ( t ) being the number of cells in population k at time t . The second equality in (3.22) is obtained byusing the normal vector − ν and the post-collision directions v (cid:48) h and v (cid:48) k in ˜˜ f hk over the set V − hk . Notice that thefollowing property holds (cid:90) V K hk ( · , · , v h ) d v h = 0 , (3.23)which ensures that the density of cells in population h will be preserved in the course of such interactions.Moreover, based on (3.20) and (3.22), we define the rate of change of the one-particle distribution function p h ( x h , t, v h ) due to population-switching interactions leading the test cell to leave population h as J hk ( x h , t, v h ) := − c | V | (cid:37) n − N k ( t ) (cid:90) V + hk (cid:90) V ν · ( v h − v k ) ˜˜ f hk ( x h , x h + ν(cid:37), t, v h , v k ) d v k d ν . (3.24)Analogously, we define the rate of change of p h ( x h , t, v h ) due to population-switching interactions leading atest cell to leave a generic population l (cid:54) = h and enter population h as J hlk ( x h , t, v h ) := c | V | (cid:37) n − N k ( t ) (cid:90) Ω l ( x k ) (cid:90) V δ ( x l − x h ) δ ( v l − v h ) × (cid:90) V + lk (cid:90) V ν · ( v l − v k ) ˜˜ f lk ( x l , x l + ν(cid:37), t, v l , v k ) d v k d ν d v l d x l , (3.25)with δ ( z − z ∗ ) being the Dirac delta distribution centred at z ∗ . Definition (3.25) ensures that the density ofcells that leave population l due to such interactions will enter population h . In fact, we have J hlk ( x h , t, v h ) = − (cid:90) Ω l ( x k ) (cid:90) V δ ( x l − x h ) δ ( v l − v h ) J lk ( x l , t, v l ) d v l d x l . In summary, the term Q hk in transport equation (3.14) is defined in terms of (3.22)-(3.25) in different possibleways depending on the cell-cell interactions that are considered.Under Assumptions 1-4 and definitions (3.22), (3.24) and (3.25), the rates of change of the one-particle distri-bution functions p I ( x I , t, v I ), p A ( x A , t, v A ) and p D ( x D , t, v D ) due to cell-cell interactions will be, respectively, Q ID ( x I , t, v I ) = (1 − ζ ) K ID ( x I , t, v I )+ ζ J ID ( x I , t, v I ) , (3.26) Q AD ( x A , t, v A ) = K AD ( x A , t, v A ) + ζ J AID ( x A , t, v A ) , (3.27) Q DI ( x D , t, v D ) = K DI ( x D , t, v D ) . (3.28)9ubstituting (3.16), (3.17) and (3.26)-(3.28) into transport equation (3.14), we obtain the following transportequations for p I ( x I , t, v I ), p A ( x A , t, v A ) and p D ( x D , t, v D ): ∂ t p I + c v I · ∇ x I p I = − ( − T I ) B [ p I ] (cid:124) (cid:123)(cid:122) (cid:125) cell motion + (1 − ζ ) K ID (cid:124) (cid:123)(cid:122) (cid:125) interactions − ζ J ID , (cid:124) (cid:123)(cid:122) (cid:125) outflow dueto activation x I ∈ R n , t ∈ R ∗ + , v I ∈ V , (3.29) ∂ t p A + c v A · ∇ x A p A = − b ( − T A )[ p A ] (cid:124) (cid:123)(cid:122) (cid:125) cell motion + K AD (cid:124) (cid:123)(cid:122) (cid:125) interactions + ζ J AID , (cid:124) (cid:123)(cid:122) (cid:125) inflow dueto activation x A ∈ R n , t ∈ R ∗ + , v A ∈ V , (3.30) ∂ t p D + c v D · ∇ x D p D = − b ( − T D )[ p D ] (cid:124) (cid:123)(cid:122) (cid:125) cell motion + K DI , (cid:124) (cid:123)(cid:122) (cid:125) interactions x D ∈ R n , t ∈ R ∗ + , v D ∈ V . (3.31)The terms on the right-hand sides of (3.29)-(3.31) represent the rate of change of the one-particle distributionsdue to the biophysical phenomena specified below each term. In this section, we derive a coupled system of balance equations for the macroscopic cell densities correspondingto the mescoscopic scale model given by transport equations (3.29)-(3.31).
Scaling
We assume the mean run time ¯ τ to be small compared to the characteristic temporal scale for thedynamics of the macroscopic cell densities, which is represented by the parameter T ∈ R ∗ + , i.e. we make theassumption ¯ τT =: ε (cid:28) . Moreover, we let X ∈ R ∗ + represent the characteristic spatial scale for the dynamics of the macroscopic celldensities and introduce the rescaled quantitiesˆ t = tT , ˆ x = x X , ˆ τ = ¯ τT , ˆ c = c TX . As similarly done in [1, 8], in order to obtain a mathematical model for the dynamics of the cells at themacroscopic scale, we consider the scaling( x , t, c, τ ) (cid:55)→ (ˆ x /ε, ˆ t/ε, ˆ c/ε γ , ˆ τ /ε µ ) , (4.1)with γ, µ ∈ R ∗ + , γ < µ > − γ . (4.2)10hroughout the rest of the paper, we will drop the carets from (4.1) and we will study two-dimensional celldynamics ( i.e. from now on we assume n = 2).Furthermore, noting that the diameter of the cells is small compared to the characteristic spatial scale forthe dynamics of the macroscopic cell densities, and considering a biological scenario where the number of cellsin the system is large and activation of CTL occurs with a small probability, we assume (cid:37) = ε ξ , N I ( t ) ≡ ε − ϑ , N D = ε − ϑ , ζ = ε κ , ξ, ϑ, κ ∈ R ∗ + . (4.3)In particular, we will be focussing on a biological scenario corresponding to the following assumptions γ := 12 , ξ − ϑ := 1 − γα − κ = − ( ξ − ϑ ) + 32 > . (4.4)Notice that ξ − ϑ < α < /
2, whereas in the case where cells undergo classical diffusion we have α = 2and, therefore, ξ − ϑ = 1 − γ = 1 /
2. Under scaling (4.1), (4.2) definitions (2.3) become ψ εi ( x i , τ i ) = (cid:16) ε µ τ ( x i ) ε µ τ ( x i ) + τ i (cid:17) α , ϕ εi ( x i , τ i ) := α ε µ τ ( x i ) α ( ε µ τ ( x i ) + τ i ) α +1 , α ∈ (1 , . (4.5)Moreover, under assumption (4.3) on (cid:37) we have˜˜ f hk ( x h , x h ± νρ, t, v h , v k ) ≡ ˜˜ f hk ( x h , x h ± ε ξ ν, t, v h , v k ) . (4.6) “Molecular chaos” assumption Considering a biological scenario where cell densities are low, we assumethe velocities of any two cells which are about to collide to be uncorrelated – i.e. we make the so-called“molecular chaos” assumption, which holds at low densities and is commonly used in kinetic theory [4, 20].Under this assumption, the two-particle distribution function ˜˜ f hk ( x h , x h ± ε ξ ν, t, v h , v k ) can be expressed asthe product of the corresponding one-particle distribution functions, that is,˜˜ f hk ( x h , x h ± ε ξ ν, t, v h , v k ) = p εh ( x h , t, v h ) p εk ( x h , t, v k ) + O ( ε ξ ) . (4.7)We draw the attention of the reader to the fact that, throughout the rest of the paper, superscript and subscript ε related to the scaling should not be confused with population indices.Under scaling (4.1), (4.2) and assumptions (4.3), using (4.6), (4.7) and assuming n = 2, the interaction termsdefined via (3.22), (3.24) and (3.25) read as K εhk ( x h , t, v h ) = 1 | V | ε ξ − ϑ − γ c (cid:90) V + hk (cid:90) V ν · ( v h − v k ) (cid:104) p εh ( x h , t, v (cid:48) h ) p εk ( x h , t, v (cid:48) k ) − p εh ( x h , t, v h ) p εk ( x h , t, v k ) (cid:105) d v k d ν , (4.8) J εhk ( x h , t, v h ) = − | V | ε ξ − ϑ − γ c (cid:90) V + hk (cid:90) V ν · ( v h − v k ) p εh ( x h , t, v h ) p εk ( x h , t, v k ) d v k d ν (4.9)and ε J hlk ( x h , t, v h ) = 1 | V | ε ξ − ϑ − γ c (cid:90) Ω l ( x k ) (cid:90) V δ ( x l − x h ) δ ( v l − v h ) × (cid:90) V + lk (cid:90) V ν · ( v l − v k ) p εl ( x l , t, v l ) p εk ( x l , t, v k ) d v k d ν d v l d x l . (4.10)11 xpansion of p εh and macroscopic cell quantities We expand the one-particle distribution function p εh in terms of its zeroth moment ρ εh ( i.e. the macroscopic cell density) and its first moment w εh ( i.e. the localmacroscopic direction of cell motion) as p εh ( x h , t, v h ) = 1 | V | (cid:16) ρ εh ( x h , t ) + ε γ v h · w εh ( x h , t ) (cid:17) + o ( ε γ ) , h ∈ { A, D, I } , (4.11)with ρ εh ( x h , t ) := (cid:90) V p εh ( x h , t, v h ) d v h , w εh ( x h , t ) := (cid:90) V v h p εh ( x h , t, v h ) d v h . (4.12) Transport equation for p εh Under scaling (4.1), (4.2) and assumptions (4.3), using (4.6), (4.7) and assuming n = 2, we rewrite transport equation (3.14) for the one-particle distribution function p h ( x h , t, v h ) as ε∂ t p εh + ε − γ c v h · ∇ x h p εh = − ( − T h )[ ε p β h h ] + Q εhk , (4.13)where Q εhk is defined in terms of K εhk , J εhk and ε J hlk as per (3.26)-(3.28), that is, Q εID ( x I , t, v I ) = (1 − ε κ ) K εID ( x I , t, v I ) + ε κ J εID ( x I , t, v I ) , (4.14) Q εAD ( x A , t, v A ) = K εAD ( x A , t, v A ) + ε κ ε J AID ( x A , t, v A ) (4.15)and Q εDI ( x D , t, v D ) = K εDI ( x D , t, v D ) . (4.16)We recall that in the case where cells move in a local search pattern ( i.e. for h = A and h = D ), β h is definedvia (2.1) and (2.2), and thus ε p β h h ( x h , t, v h ) is given as in (3.16). On the other hand, in the case where cellsmove in a non-local search pattern ( i.e. for h = I ), β h is defined via (2.1) and (2.3), and thus ε p β h h ( x h , t, v h ) isgiven by (3.17) with B ε [ p εh ]( x h , t, v h ) = (cid:90) t B ε ( x h , t − s ) p εh ( x h − ( c v h + b )( t − s ) , s ) d s . As before, B ε is defined through its Laplace transform in time ˆ B ε and, in particular, under assumptions (4.2),we make the approximationˆ B ε ( x h , ελ + ε µ b + ε − γ c v h · ∇ x h ) (cid:39) ˆ B ε ( x h , ε − γ c v h · ∇ x h ) . Using the properties of the Laplace transform of a convolution, we write (cid:90) t B ε ( x h , t − s ) p εh ( x h − ( c v h + b )( t − s ) , s, v h ) d s (cid:39) ˆ B ε ( x h , ε − γ c v h · ∇ x h ) p εh ( x h , t, v h ) , with ˆ B ε ( x h , ε − γ c v h · ∇ x h ) = ˆ ϕ εh ( x h , ε − γ c v h · ∇ x h )ˆ ψ εh ( x h , ε − γ c v h · ∇ x h ) . (4.17)Substituting the expressions for ˆ ϕ εh and ˆ ψ εh into (4.17), calculations similar to those carried out in [8, 9] allowone to show that ˆ B ε ( x h , ε − γ c v h · ∇ x h ) = α − d ε − ε − γ − α c v h · ∇ x h − d α − ε ε (1 − γ )( α − ( c v h · ∇ x h ) α − ( α − Γ( − α + 1) + O ( d α − ε λ α ) . (4.18)In (4.18), d ε ( x h ) := τ ( x h ) ε µ , where τ ( x h ) is defined via (2.3), and Γ( · ) denotes the gamma function.12 ransport equations for ρ εI , ρ εA and ρ εD Integrating both sides of transport equation (4.13) with respect to v h over the set V and using the fact that the turning operator T h satisfies (2.6), we find that the macroscopiccell density ρ εh ( x h , t ) given by (4.12) satisfies the following transport equation ∂ t ρ εh + 2 c ∇ x h · w εh = ε − (cid:90) V Q εhk d v h , x h ∈ R n , t ∈ R ∗ + . (4.19)Moreover, substituting the expressions for p εh ( x h , t, v h ) and p εk ( x h , t, v k ) given by (4.11) into the definitionsof K εhk , J εhk and ε J lhk given by (4.8)-(4.10) we find (cid:90) V K εhk ( x h , t, v h ) d v h = 0 (4.20)and, neglecting higher order terms, (cid:90) V J εhk ( x h , t, v h ) d v h = − ε ξ − ϑ − γ c M ρ εh ρ εk , (cid:90) V ε J hlk ( x h , t, v h ) d v h = ε ξ − ϑ − γ c M ρ εl ρ εk , (4.21)where M is defined as M := 1 | V | (cid:90) V (cid:90) V (cid:90) V + hk ν · ( v h − v k ) d ν d v h d v k , h, k ∈ { A, D, I } , h (cid:54) = k . (4.22)Notice that relation (4.20) is obtained using property (3.23).In conclusion, using (4.20) and (4.21) along with (4.14)-(4.16) and the scaling assumption on κ given by (4.4),from transport equation (4.19) we obtain the following balance equations for the macroscopic cell densities ρ εI ( x I , t ), ρ εA ( x A , t ) and ρ εD ( x D , t ) ∂ t ρ εI + 2 c ∇ x I · w εI = − c M ρ εI ρ εD , x I ∈ R , t ∈ R ∗ + , (4.23) ∂ t ρ εA + 2 c ∇ x A · w εA = c M ρ εI ρ εD , x A ∈ R , t ∈ R ∗ + , (4.24) ∂ t ρ εD + 2 c ∇ x D · w εD = 0 , x D ∈ R , t ∈ R ∗ + . (4.25)The terms on the right-hand sides of (4.23) and (4.24) model the rate of change of the cell densities ρ εI and ρ εA ,respectively, due to the fact that immune activation via interactions with cells of population D leads cells toleave population I and enter population A . Transport equations for w εI , w εA and w εD Multiplying both sides of transport equation (4.13) by v h andthen integrating both sides of the resulting equation with respect to v h over the set V, we find that the localmacroscopic direction of cell motion w εh ( x h , t ) given by (4.12) satisfies the following equation ε γ ∂ t w εh + ε − γ c ∇ x h (cid:90) V v h ⊗ v h p εh d v h = − (cid:90) V v h ( − T h )[ ε p β h h ] d v h + (cid:90) V v h Q εhk d v h . (4.26)In the case where β h is defined via (2.1) and (2.2), using (4.11), (3.16) and the properties of the turningoperator T h established by Lemma 1 in Appendix D, we find that the first term on the right-hand side of (4.26)is given by (cid:90) V v h ( − T h )[ ε p β h h ] d v h = 2 ε γ | V | b ( ι − w εh . (4.27)13ere, ι ( x h , t ) is the first non-zero eigenvalue of the turning operator T h , which is given by (D.1). On the otherhand, when β h is defined via (2.1) and (2.3), using (4.11), (3.17) and the properties of the turning operator T h established by Lemma 1 in Appendix D, it was proved in [9] that the following approximate expression for thefirst term on the right-hand side of (4.26) holds (cid:90) V v h ( − T h )[ ε p β h h ] d v h = ε − γα − (cid:16) g α ∇ α − ρ εh + 2( α − τ | V | ( ι − w εh (cid:17) + l.o.t. , (4.28)where g α ( x h , t ) := − πτ α − (1 − α ) c α − sin( πα )Γ( α ) (cid:16) ι − | V || V | (cid:17) for Γ( − α + 1) = π sin( πα )Γ( α ) . (4.29)Notice that g α ( · , · ) > πα ) < α ∈ (1 , · ) (cid:48) : V (cid:55)→ V is a bijection and v (cid:48) h · ν = − v h · ν , whence ν · ( v h − v k ) = − ν · ( v (cid:48) h − v (cid:48) k ), we find (cid:90) V v h K εhk d v h == 1 | V | ε ξ − ϑ − γ c (cid:16) (cid:90) V (cid:90) V (cid:90) V + hk v h p εh ( x h , t, v (cid:48) h ) p εk ( x h , t, v (cid:48) k ) ν · ( v h − v k ) d ν d v k d v h − (cid:90) V (cid:90) V (cid:90) V + hk v h p εh ( x h , t, v h ) p εk ( x h , t, v k ) ν · ( v h − v k ) d ν d v k d v h (cid:17) = − | V | ε ξ − ϑ − γ c (cid:90) V (cid:90) V (cid:90) V + hk ( v (cid:48) h ) (cid:48) p εh ( x h , t, v (cid:48) h ) p εk ( x h , t, v (cid:48) k ) ν · ( v (cid:48) h − v (cid:48) k ) d ν d v (cid:48) h d v (cid:48) k − (cid:90) V (cid:90) V (cid:90) V + hk v h p εh ( x h , t, v h ) p εk ( x h , t, v k ) ν · ( v h − v k ) d ν d v h d v k = 1 | V | ε ξ − ϑ − γ c (cid:90) V (cid:90) V (cid:90) V + hk ( v (cid:48) h − v h ) p εh ( x h , t, v h ) p εk ( x h , t, v k ) ν · ( v h − v k ) d ν d v h d v k = − | V | ε ξ − ϑ − γ c (cid:90) V (cid:90) V | v h − v k | v h p εh ( x h , t, v h ) p εk ( x h , t, v k ) d v h d v k . (4.30)The last equality in (4.30) is obtained using the fact that v (cid:48) h − v h = − v h · ν ) ν . Substituting the expressionsfor p εh ( x h , t, v h ) and p εk ( x h , t, v k ) given by (4.11) into (4.30) and into definitions (4.9) and (4.10), neglectinghigher order terms we find (cid:90) V v h K εhk d v h = − ε ξ − ϑ c
83 1 | V | q h ρ εk w εh , (4.31)and (cid:90) V v h J εhk d v h = − ε ξ − ϑ c M ρ εk w εh , (cid:90) V v h ε J hlk d v h = ε ξ − ϑ c M ρ εk w εh , (4.32)where M is given by (4.22) and q h is defined as q h := (cid:90) V | v h − v k | d v , h ∈ { A, D, I } . (4.33)Finally, substituting the expression for p εh ( x h , t, v h ) into the second term on the left-hand side of (4.26) andneglecting higher order terms yields c ∇ x h (cid:90) V v h ⊗ v h p εh d v h = C h ∇ x h ρ εh , (4.34)14ith C h being defined as C h := c | V | (cid:90) V v h ⊗ v h d v h , h ∈ { A, D, I } . (4.35)In conclusion, using (3.16) and (3.17), (4.31), (4.32) and (4.34) along with (4.14)-(4.16), from equation (4.26)we obtain the following equations for the local macroscopic directions of cell motion w εI ( x I , t ), w εA ( x A , t ) and w εD ( x D , t ): ε γ ∂ t w εI + ε − γ C I ∇ x I ρ εI = ε − γα − (cid:18) g α ∇ α − ρ εI + 2( α − τ | V | ( ι − w εI (cid:19) − ε ξ − ϑ c (cid:18) (1 − ε κ ) 83 1 | V | q I ρ εD w εI + 2 ε κ M ρ εD w εI (cid:19) , x I ∈ R , t ∈ R ∗ + , (4.36) ε γ ∂ t w εA + ε − γ C A ∇ x A ρ εA = − ε γ | V | b ( ι − w εA − ε ξ − ϑ c (cid:18)
83 1 | V | q A ρ εD w εA − ε κ M ρ εD w εI (cid:19) , x A ∈ R , t ∈ R ∗ + , (4.37)and ε γ +1 ∂ t w εD + ε − γ C D ∇ x D ρ εD = − ε γ | V | b ( ι − w εD − ε ξ − ϑ c
83 1 | V | q D ρ εI w εD , x D ∈ R , t ∈ R ∗ + . (4.38) Macroscopic scale model
Noting that 1 − γα − < − γ since α ∈ (1 ,
2) and using assumptions (4.4),letting ε → w I ( x , t ), w A ( x , t ) and w D ( x , t ) of the asymptotic expansions for the local macroscopic directions of cell motion w Iε ( x , t ), w Aε ( x , t ) and w Dε ( x , t ) w I = − g α σ α + σ q I ρ D ∇ α − x ρ I , w A = − C A σ b − σ q A ρ D ∇ x ρ A , w D = − C D σ b − σ q D ρ I ∇ x ρ D (4.39)with σ α ( x , t ) := 2( α − ι ( x , t ) − τ ( x ) | V | , σ b ( x , t ) := 2 b | V | ( ι ( x , t ) − ,σ q h := − c | V | q h , for h ∈ { I, A, D } . (4.40)Furthermore, under assumptions (4.4), letting ε → ρ I ( x , t ), ρ A ( x , t ) and ρ D ( x , t ) of the asymptoticexpansions for the macroscopic cell densities ρ εI ( x , t ), ρ εA ( x , t ) and ρ εD ( x , t ) ∂ t ρ I − c ∇ x · (cid:16) g α σ α + σ q I ρ D ∇ α − x ρ I (cid:17) = − c M ρ I ρ D , α ∈ (1 , , x ∈ R , t ∈ R ∗ + , (4.41) ∂ t ρ A − c ∇ x · (cid:16) C A σ b − σ q A ρ D ∇ x ρ A (cid:17) = c M ρ I ρ D , x ∈ R , t ∈ R ∗ + , (4.42) ∂ t ρ D − c ∇ x · (cid:16) C D σ b − σ q D ρ I ∇ x ρ D (cid:17) = 0 , x ∈ R , t ∈ R ∗ + , (4.43)15here M is defined via (4.22) and g α is defined according to (4.29). Notice that the dependence of the diffusioncoefficients in (4.41)-(4.43) on the cell densities follows from conservative interactions between cells of differentpopulations, while population-switching interactions do not affect the diffusion coefficients. The modelling approach for the switch between cell migration modes presented here could be generalised byincluding additional cellular phenomena involved in the immune response to cancer, and considering otheraspects of immune cell movement as well.With reference to the mathematical modelling of the immune response to cancer, a natural generalisationwould be to include a population of cancer cells and allow activated CTLs to induce death in cancer cellsvia binary interactions. Moreover, the recognition phase of the adaptive immune response to cancer could bemodelled by splitting the population of DCs into a subpopulation of cells with no tumour antigens on theirsurface and a subpopulation of cells presenting some antigen – which would move in a non-local and in a morelocalised search pattern, respectively [7] – and letting DCs switch from one subpopulation to the other viabinary interactions with cancer cells [21, 22]. The strategy we have used here to model population-switchinginteractions may prove useful to the development of both generalisations of our modelling approach.In regard to the mathematical modelling of other aspects of immune cell movement, our modelling approachcould be extended to represent other switches in T cell migration patterns observed in the immune responseto different pathogens, which are driven by possible chemotactic cues and by the conditions of the surroundingmicroenvironment [13]. Moreover, further generalisations of the modelling approach could be developed inrelation to experimental results indicating that T cells can also undergo subdiffusive [25] and fully ballistic [24]migration.In general, it would be interesting to apply the modelling approach presented in this paper and its possibledevelopments to other biological and ecological contexts whereby switch from non-local to localised migrationpatterns has been reported [2, 5, 6, 11, 16]. A recent work in this direction is [19], where a model for the switchbetween L´evy and Brownian movement determined by internal chemical pathways in bacteria was considered.
A Derivation of transport equations (3.11) - (3.13) Using the method presented in [8], we show how to derive a transport equation for the two-particle distributionfunction ˜˜ f hk ( x h , x k , t, v h , v k ) starting from transport equation (3.2) for the two-particle distribution function f hk ( x h , x k , t, v h , v k , τ h , τ k ).We first introduce the notation˜ f τ h ( x h , x k , t, v h , v k , τ h ) := (cid:90) t f hk ( x h , x k , t, v h , v k , τ h , τ k ) d τ k , and ˜ f τ k ( x h , x k , t, v h , v k , τ k ) := (cid:90) t f hk ( x h , x k , t, v h , v k , τ h , τ k ) d τ h , and then note that, when β h and β k are given by (2.1) with ψ h and ψ k defined via (2.2) or (2.3), the solutionsof (3.2) subject to the initial and boundary conditions considered here are such that ˜ f τ h decays monotonicallyas τ h increases, and ˜ f τ k exhibits an analogous behaviour. Hence, integrating (3.2) with respect to ( τ h , τ k ) over(0 , t ) with t large enough so that ˜ f τ h ( x h , x k , t, v h , v k , τ h = t ) is negligible compared to ˜ f τ h ( x h , x k , t, v h , v k )16nd ˜ f τ k ( x h , x k , t, v h , v k , τ k = t ) is negligible compared to ˜ f τ k ( x h , x k , t, v h , v k ), with˜ f τ h ( x h , x k , t, v h , v k ) =: ˜ f τ h ( x h , x k , t, v h , v k , τ h = 0)and ˜ f τ k ( x h , x k , t, v h , v k ) =: ˜ f τ k ( x h , x k , t, v h , v k , τ k = 0) , we obtain the following transport equation for ˜˜ f hk ( x h , x k , t, v h , v k )( ∂ t + c v h · ∇ x h + c v k · ∇ x k ) ˜˜ f hk = − (cid:90) t (cid:90) t ( β h + β k ) f hk d τ h d τ k + ˜ f τ h + ˜ f τ k , which can be rewritten as( ∂ t + c v h · ∇ x h + c v k · ∇ x k ) ˜˜ f hk = − ˜˜ f β h hk − ˜˜ f β k hk + ˜ f τ h + ˜ f τ k , (A.1)with ˜˜ f β h hk ( x h , x k , t, v h , v k ) and ˜˜ f β k hk ( x h , x k , t, v h , v k ) given by (3.6).When cell movement at the microscopic scale obeys the rules presented in Section 2, we have˜ f τ h = T h [ ˜˜ f β h hk ] and ˜ f τ k = T k [ ˜˜ f β k hk ] , (A.2)with the turning operators T h and T k being defined via (2.5). The first two terms on the right-hand side of (A.1)describe the density of cells that stop with rates β h , β k . The initial conditions at τ h = 0 and τ k = 0 ( i.e. atthe beginning of a new run phase) given by (A.2) describes how the cells will resume their motion in a newdirection dictated by the turning operators T h and T k , respectively.Substituting the expressions for ˜ f τ h and ˜ f τ k given by (A.2) into transport equation (A.1) yields( ∂ t + c v h · ∇ x h + c v k · ∇ x k ) ˜˜ f hk = − ( − T h )[ ˜˜ f β h hk ] − ( − T k )[ ˜˜ f β k hk ] , ( x h , x k ) ∈ Ω , t ∈ R + , ( v h , v k ) ∈ V . (A.3) Remark A.1.
Since we consider transport equation (3.2) complemented with a smooth, compactly supportedinitial condition, the initial condition for transport equation (A.3) will be a smooth, compactly supported functionas well. Therefore, the two-particle distribution function ˜˜ f hk ( x h , x k , t, v h , v k ) will have compact support on Ω × V for all t ∈ R ∗ + . B Derivation of the equation for the one-particle distribution
Transport equation (3.14) for the one-particle distribution function p h ( x h , t, v h ) can derived from transportequation (A.3) for the two-particle distribution ˜˜ f hk ( x h , x k , t, v h , v k ) in six steps as previously done in [8]. (I) We integrate transport equation (A.3) with respect to ( x k , v k ) over the set Ω k ( x h ) × V and multiplyboth sides of the resulting equation by | V | − to obtain | V | − (cid:90) Ω k ( x h ) (cid:90) V ( ∂ t + c v h · ∇ x h + c v k · ∇ x k ) ˜˜ f hk d v k d x k = − | V | − (cid:90) Ω k ( x h ) (cid:90) V ( − T h )[ ˜˜ f β h hk ] d v k d x k − | V | − (cid:90) Ω k ( x h ) (cid:90) V ( − T k )[ ˜˜ f β k hk ] d v k d x k . (B.1)17 II)
Using the fact that p h is given by (3.5) and integrals with respect to x k and v k commute, we rewrite thefirst term on the left-hand side of (B.1) as | V | − ∂ t (cid:90) Ω k ( x h ) (cid:90) V ˜˜ f hk d v k d x k = ∂ t p h . (III) Using Reynold’s transport theorem in the variable x h , we rewrite the second term on the left-hand sideof (B.1) as | V | − c (cid:90) Ω k ( x h ) (cid:90) V ( v h · ∇ x h ) ˜˜ f hk d v k d x k = | V | − c v h · ∇ x h p h − | V | − c (cid:90) ∂ B (cid:37) ( x h ) (cid:90) V ( v h · ν ) ˜˜ f hk d v k d σ . Here, ν is the unit normal to ∂ Ω k ( x h ) that points outward from Ω k ( x h ) and inward to B (cid:37) ( x i ), and d σ denotesthe surface element. (IV) Since ˜˜ f hk has compact support on Ω × V ( vid. Remark A.1), we use the divergence theorem andrewrite the third term on the left-hand side of (B.1) as | V | − c (cid:90) Ω k ( x h ) (cid:90) V ( v k · ∇ x k ) ˜˜ f hk d v k d x k = | V | − c (cid:90) ∂ B (cid:37) ( x h ) (cid:90) V ( v k · ν ) ˜˜ f hk d v k d σ . (V) Changing order of integration, we rewrite the first term on the right-hand side of (B.1) as −| V | − (cid:90) Ω k ( x h ) (cid:90) V ( − T h )[ ˜˜ f β h hk ] d v k d x k = − ( − T h )[ p β h h ] , with p β h h ( x h , t, v h ) given by (3.7). (VI) Since T k satisfies (2.6), the second term on the right-hand side of (B.1) is identically zero.Taken together, the results obtained in Steps (I) - (VI) allow one to conclude that the one-particle distributionfunction p h ( x h , t, v h ) satisfies the following transport equation ∂ t p h + c v h · ∇ x h p h = − ( − T h )[ p β h h ] + Q hk , x h ∈ R n , t ∈ R + , v h ∈ V , (B.2)with the weighted one-particle distribution function p β h h ( x h , t, v h ) being given by (3.7) and the term Q hk ( x h , t, v h )being defined according to (3.15). C Derivation of the non-local trajectory term
Here we show how to derive the non-local trajectory term (3.17). In the case where β h is defined via (2.1)and (2.3) ( i.e. for h = I ) and β k is defined via (2.1) and (2.2) ( i.e. for k = D ), applying the method ofcharacteristics to (3.2) and using the fact that ψ k ( · , τ k ) ψ k ( · , τ k − τ h ) = e − bτ h one finds [8] f hk = f hk ( x h − c v h τ h , x k − c v k τ h , t − τ h , v h , v k , τ h = 0 , τ k − τ h ) ψ h ( x h , τ h ) e − bτ h . (C.1)18ntroducing the notation¯ f hk := (cid:90) Ω k ( x h ) (cid:90) V (cid:90) t f hk ( · , x k − c v k τ h , · , · , v k , τ h = 0 , τ k − τ h ) d τ k d v k d x k and substituting (C.1) into (3.7) gives p β h h ( x h , t, v h ) = 1 | V | (cid:90) t ϕ h ( x h , τ h ) ψ h ( x h , τ h ) (cid:90) Ω k ( x h ) (cid:90) V (cid:90) t f hk d τ k d v k d x k d τ h = 1 | V | (cid:90) t ϕ h ( x h , τ h ) e − bτ h ¯ f hk ( x h − c v h τ h , t − τ h , v h ) d τ h = 1 | V | (cid:90) t ϕ h ( x h , t − s ) e − ( t − s )( b + c v h ·∇ x h ) ¯ f hk ( x h , s, v h ) d s . (C.2)The last equality in (C.2) is obtained using the change of variables s = t − τ h along with the following Taylorexpansion e − ( t − s ) c v ·∇ f ( x ) = ∞ (cid:88) m =0 ( − ( t − s ) c v · ∇ ) m m ! f ( x )= ∞ (cid:88) m =0 m ! ( − ( t − s ) c v ) m ∇ m f ( x ) = f ( x − ( t − s ) c v ) . Hence, the Laplace transform in time of p β h h ( x h , t, v h ) isˆ p β h h ( x h , λ, v h ) = 1 | V | ˆ ϕ h ( x h , λ + b + c v h · ∇ x h ) ˆ¯ f hk ( x h , λ, v h ) . (C.3)Here, λ is the Laplace variable, and ˆ ϕ h and ˆ¯ f hk are the Laplace transforms in time of the functions ϕ h and ¯ f hk .Moreover, substituting (C.1) into (3.5) and computing the Laplace transform in time yieldsˆ p h ( x h , λ, v h ) = 1 | V | ˆ ψ h ( x h , λ + b + c v h · ∇ x h ) ˆ¯ f hk ( x h , λ, v h ) , with ˆ ψ h being the Laplace transform of the function ψ h . The latter equation givesˆ¯ f hk ( x h , λ, v h ) = | V | ˆ p h ( x h , λ, v h )ˆ ψ h ( x h , λ + b + c v h · ∇ x h ) . Substituting such an expression for ˆ¯ f hk into (C.3) ones sees that (C.2) can be written as p β h h ( x h , t, v h ) = B [ p h ]( x h , t, v h )with the integral operator B being defined according to (3.18). D Properties of the turning operator T h Let v ≡ ( v , v , . . . , v n − ) ∈ V. Moreover, denote by | V | the surface area of the unit n -sphere V, i.e. | V | = π n / Γ (cid:0) n (cid:1) , for n even ,π n / Γ (cid:0) n + 1 (cid:1) , for n odd , e = (1 , , . . . , ∈ R n the vector with all components equal to 1. Some useful properties of the spectrumof the turning operator T h defined via (2.5), where the turning kernel (cid:96) satisfies (2.4), are established by thefollowing lemma [1]. Lemma 1.
If the turning kernel (cid:96) ( · , · , | ¯ v − v | ) is continuous, then T h is a symmetric compact operator. Inparticular, there exists an orthonormal basis of L (V) consisting of eigenfunctions of T h and: φ ( v ) = 1 | V | is an eigenfunction of T h with eigenvalue ι = 1 ; φ i ( v ) = n v i | V | for i = 0 , . . . , n − are eigenfunctions of T h with eigenvalue ι ( · , · ) = (cid:90) V (cid:96) ( · , · , | v − e | ) v d v < . (D.1) References [1] Wolgang Alt. Biased random walk models for chemotaxis and related diffusion approximations.
Journalof Mathematical Biology , 9(2):147–177, 1980.[2] Frederic Bartumeus, Francesc Peters, Salvador Pueyo, Celia Marras´e, and Jordi Catalan. Helical l´evy walks:adjusting searching statistics to resource availability in microzooplankton.
Proceedings of the NationalAcademy of Sciences , 100(22):12771–12775, 2003.[3] Alexandre Boissonnas, Luc Fetler, Ingrid S Zeelenberg, St´ephanie Hugues, and Sebastian Amigorena. Invivo imaging of cytotoxic t cell infiltration and elimination of a solid tumor.
The Journal of ExperimentalMedicine , 204(2):345–356, 2007.[4] Carlo Cercignani, Reinhard Illner, and Mario Pulvirenti.
The mathematical theory of dilute gases , volume106. Springer Science & Business Media, 2013.[5] Monique de Jager, Frederic Bartumeus, Andrea K¨olzsch, Franz J Weissing, Geerten M Hengeveld, Bart ANolet, Peter MJ Herman, and Johan van de Koppel. How superdiffusion gets arrested: ecological encountersexplain shift from l´evy to brownian movement.
Proceedings of the Royal Society B: Biological Sciences ,281(1774):20132605, 2014.[6] HJ De Knegt, GM Hengeveld, F Van Langevelde, WF De Boer, and KP Kirkman. Patch density determinesmovement patterns and foraging efficiency of large herbivores.
Behavioral Ecology , 18(6):1065–1072, 2007.[7] John J Engelhardt, Bijan Boldajipour, Peter Beemiller, Priya Pandurangi, Caitlin Sorensen, Zena Werb,Mikala Egeblad, and Matthew F Krummel. Marginating dendritic cells of the tumor microenvironmentcross-present tumor antigens and stably engage tumor-specific t cells.
Cancer Cell , 21(3):402–417, 2012.[8] Gissell Estrada-Rodriguez and Heiko Gimperlein. Interacting particles with levy strategies: limits oftransport equations for swarm robotic systems.
SIAM Journal on Applied Mathematics , 80(1):476–498,2020.[9] Gissell Estrada-Rodriguez, Heiko Gimperlein, and Kevin J Painter. Fractional patlak–keller–segel equationsfor chemotactic superdiffusion.
SIAM Journal on Applied Mathematics , 78(2):1155–1173, 2018.2010] Benjamin Franz, Jake P Taylor-King, Christian Yates, and Radek Erban. Hard-sphere interactions invelocity-jump models.
Physical Review E , 94(1):012129, 2016.[11] Nicolas E Humphries, Nuno Queiroz, Jennifer RM Dyer, Nicolas G Pade, Michael K Musyl, Kurt MSchaefer, Daniel W Fuller, Juerg M Brunnschweiler, Thomas K Doyle, Jonathan DR Houghton, et al.Environmental context explains l´evy and brownian movement patterns of marine predators.
Nature ,465(7301):1066–1069, 2010.[12] Earle H Kennard et al.
Kinetic theory of gases , volume 287. McGraw-hill New York, 1938.[13] Matthew F Krummel, Frederic Bartumeus, and Audrey G´erard. T cell migration, search strategies andmechanisms.
Nature Reviews Immunology , 16(3):193, 2016.[14] Fiona R Macfarlane, Mark AJ Chaplain, and Tommaso Lorenzi. A stochastic individual-based modelto explore the role of spatial interactions and antigen recognition in the immune response against solidtumours.
Journal of Theoretical Biology , 480:43–55, 2019.[15] Fiona R Macfarlane, Tommaso Lorenzi, and Mark AJ Chaplain. Modelling the immune response to cancer:an individual-based approach accounting for the difference in movement between inactive and activated tcells.
Bulletin of Mathematical Biology , 80(6):1539–1562, 2018.[16] Bart A Nolet and Wolf M Mooij. Search paths of swans foraging on spatially autocorrelated tubers.
Journalof Animal Ecology , pages 451–462, 2002.[17] Hans G Othmer, Steven R Dunbar, and Wolfgang Alt. Models of dispersal in biological systems.
Journalof Mathematical Biology , 26(3):263–298, 1988.[18] Hans G Othmer and Thomas Hillen. The diffusion limit of transport equations derived from velocity-jumpprocesses.
SIAM Journal on Applied Mathematics , 61(3):751–775, 2000.[19] Benoˆıt Perthame, Weiran Sun, and Min Tang. The fractional diffusion limit of a kinetic model withbiochemical pathway.
Zeitschrift f¨ur angewandte Mathematik und Physik , 69(3):67, 2018.[20] C´edric Villani. A review of mathematical topics in collisional kinetic theory.
Handbook of mathematicalfluid dynamics , 1(71-305):3–8, 2002.[21] Alex D Waldman, Jill M Fritz, and Michael J Lenardo. A guide to cancer immunotherapy: from t cellbasic science to clinical practice.
Nature Reviews Immunology , pages 1–18, 2020.[22] Stefanie K Wculek, Francisco J Cueto, Adriana M Mujal, Ignacio Melero, Matthew F Krummel, and DavidSancho. Dendritic cells in cancer immunology and immunotherapy.
Nature Reviews Immunology , pages1–18, 2019.[23] E John Wherry and Makoto Kurachi. Molecular and cellular insights into t cell exhaustion.
Nature ReviewsImmunology , 15(8):486–499, 2015.[24] Colleen M Witt, Subhadip Raychaudhuri, Brian Schaefer, Arup K Chakraborty, and Ellen A Robey.Directed migration of positively selected thymocytes visualized in real time.
PLOS Biolology , 3(6):e160,2005.[25] Tim Worbs, Thorsten R Mempel, Jasmin B¨olter, Ulrich H von Andrian, and Reinhold F¨orster. Ccr7ligands stimulate the intranodal motility of t lymphocytes in vivo.