Mass-preserving approximation of a chemotaxis multi-domain transmission model for microfluidic chips

The present work was inspired by the recent developments in laboratory experiments made on chip, where culturing of multiple cell species was possible. The model is based on coupled reaction-diffusion-transport equations with chemotaxis, and takes into account the interactions among cell populations and the possibility of drug administration for drug testing effects. Our effort was devoted to the development of a simulation tool that is able to reproduce the chemotactic movement and the interactions between different cell species (immune and cancer cells) living in microfluidic chip environment. The main issues faced in this work are the introduction of mass-preserving and positivity-preserving conditions involving the balancing of incoming and outgoing fluxes passing through interfaces between 2D and 1D domains of the chip and the development of mass-preserving and positivity preserving numerical conditions at the external boundaries and at the interfaces between 2D and 1D domains.


Introduction
The aim of the present work is to study the modelling and numerics of a chemotaxis-reaction-diffusion mathematical model describing the qualitative behavior of different cell species living in a confined environment. This work was inspired by laboratory experiments made on microfluidic chip [38], where some populations cohexist and interact. In recent years, indeed, there was the development of a new approach to biological studies aimed at reconstructing organs and complex biological processes on-chip [7]. The fundamental idea is that the comprehension of the sophisticated physiology of organisms, based on the complex behavior and interaction of cell populations, tissues and organs, needs interdisciplinary contributions, from biology to mathematics.
Motivated by laboratory setting of the experiment in microfluidic chips [7,32,38], we introduce a model describing the interactions between two cells populations, namely immune and cancer cells. The microfluidic chip is represented as a network of channels connecting two boxes (the microfluidic chambers), see Fig. 1 and a schematic picture of the experiment in Fig. 2. The mathematical model, proposed in section 3, is a reaction-diffusion system with chemotaxis and it describes birth and death processes, migration of immune cells driven by chemical signals produced by tumor cells, interaction between different cell species. From the mathematical point of view, we follow the framework of the classical macroscopic models of chemotaxis, where the evolution of the density of cells is described by a parabolic equation and the concentration of a chemoattractant can be given by a parabolic or elliptic equation, depending on the different regimes to be described and on authors' choices. The choice of a continuous model to reproduce an experiment in a confined environment, with a relatively small number of cells, is motivated by the fact that we aim at developing a simulation tool which is able to describe the phenomena of immunosorveillance of cancer in tissues, where billions of cells are present. For this reason a macroscopic model is more suitable respect a particle model. In the chambers we consider a 2D doubly-parabolic model which is a modification of the Keller-Segel model [23] to take into account the presence of two populations both producing chemical signal which are interacting each other. We remark that we consider only the 2D case since the experimental data do not take into account the height of the chip. Clearly, in principle, our framework could be easily extended to the third dimension. For the 1D microchannels connecting the 2D chambers we choose two different approaches: we can assign a 1D version of doubly-parabolic model used in the chambers; otherwise, we can assign a model derived from 1D-GA model [19], being characterized by the more realistic feature that the speed of propagation of cells in the channels is finite, which seems the dominant property at this scale. On the other hand, other models based on hyperbolic/kinetic equations for the evolution of the density of individuals can be assigned, characterized by a finite speed of propagation [16,31,15,13,12].
1.1. Original contribution of the present paper. From the mathematical and numerical viewpoint, here we deal with a challenging issue arising in chemotaxis modelling of cell interaction. The problem involves doubly-parabolic models in 2D domains (microfluifidic chambers) that are connected with 1D domains represented by channels, where either a doubly-parabolic or a hyperbolic-parabolic model can be assigned. The classical doubly-parabolic Keller-Segel (KS) model [23] of chemotaxis reads as: with u the density of individuals in the considered medium, ν the diffusion rate of the organism according to Fick's Law and φ the density of chemoattractant. The positive constant D is the diffusion coefficient of the chemoattractant; the positive coefficients a and b, are respectively its production and degradation rates, and χ is the chemotactic sensitivity, depending on the density of the considered quantities.
In the 2D domains given by the microfluidic chambers we apply a reaction-diffusion chemotaxis KS-like model inspired by (1) and described in 3.1.
In the 1D microfluidic channels, we use the one-dimensional version of the KSlike model used in the chambers, but we also studied the behavior of individuals when a hyperbolic-parabolic model, characterized by finite speed of propagation is assigned. Such hyperbolic-parabolic model, described in 3.1, is inspired by the Greeberg-Alt (GA) model, arising as a simple model for chemotaxis on a line: Note that here v is the averaged flux. Let us underline that the flux v in model (2) corresponds to v = −λ 2 ∇u + χ(u, φ)∇φ for the KS system. This system was analytically studied on the whole line and on bounded intervals in [20], while an effective numerical approximation, the Asymptotic High Order (AHO) scheme, was introduced in [28], see also [17,18] and extended on networks with general boundary conditions in [5] and [6].
Here, in the numerical treatment for the computation of numerical solutions one has to take care of what happens at the inner boundaries with the switching from 2D-doubly-parabolic models and 1D-doubly-parabolic or 1D-hyperbolic-parabolic ones.
Since we aim at reproducing the numerical solutions of such models, we need to deal with a multi-domain problem given by the passage from a 2D domain represented by the chambers of the chip to 1D domains given by the channels. For this reason we need to develop ad hoc transmission conditions to ensure the mass conservation at the 2D-1D interfaces. From the numerical viewpoint, here we consider numerical boundary conditions including in the stencil a ghost cell value taken from the neighbouring domain, as we will show in the numerical Section 4. The approximation of doubly-parabolic chemotaxis models for the 1D-KS model (1) on networks was already considered in [4]. However, in that case the transmission conditions were between 1D-1D interfaces and on each arc of the network the same fully-parabolic model was considered. We also underline that in such work, transmission conditions require to impose the continuity of the density of both cells u and chemoattractant φ, while we only impose the continuity of the fluxes, which seems to be more realistic when dealing with flux of individuals or molecules. For the numerical approximation of the GA system (2), we refer to our previous papers [28] for a single line, where the numerical treatment of the hyperbolic part of the system was based on the AHO scheme with the development of mass-preserving numerical scheme at outer boundaries, while the parabolic part was approximated by finite difference and Crank-Nicolson scheme. In [5] and [6] the GA system was solved on networks, thus making necessary to develop mass-preserving transmission conditions at inner nodes and suitable boundary conditions at outer nodes. However, the study of transmission conditions only involved the mass exchange between 1D-1D interfaces; moreover, on each arc of the network the same model was considered. Furthermore, the second order numerical approximation of the boundary conditions developed in such papers did not ensure the posivity preserving property in case of obscillating functions.
The numerical approximation of permeability Kedem-Katchalsky [22] conditions describing the conservation of the flux through a node was already considered in [34], but we underline that in the mentioned paper the study was done for the approximation with finite elements methods of linear problems. For reaction-diffusion problems the approximation of permeability conditions was studied in [36] for finite difference schemes and in [8] for discontinuous Galerkin methods. The numerical treatment of permeability conditions for chemotaxis problems was presented for the first time in [11] for the 1D parabolic-parabolic interface, and a finite difference approximation was developed without taking into consideration the mass-preservation nor the positivity-preservation properties at inner nodes. Therefore, to our knowledge, the present paper is the first numerical work where this new technique of switching size of the domains and type of equations (parabolic vs hyperbolic approach) is introduced, in order to develop mass-preserving and positivity preserving schemes.
1.2. Main contents and plan of the paper. In the present paper, a positivitypreserving and mass-preserving numerical discretization of Neumann boundary conditions at the corners and at the bottom and top boundaries of the 2D domain for 2D-doubly parabolic reaction-diffusion problem are presented. Moreover, a positivity-preserving and mass-preserving numerical scheme at the inner nodes of the network connecting the 2D chambers with the 1D channels (where the 1Ddoubly-parabolic or 1D-hyperbolic-parabolic problem can be assigned) is developed. To summarize the main contents of the present work, the mathematical issues faced in this study are indentified into two aspects: • the study of the behavior of two different modelling of the dynamics in the channels: the parabolic model describing the dynamics inside the chambers was coupled both with KS-like and GA-like models; • the numerical approximation of equations defined in a heterogeneous domain, characterized by the switch from 2D domains, represented by microfluidic left and right chambers, to 1D domains, given by the channels connecting them. Then, the numerical questions arising in the mentioned issues and here addressed are: • the study of positivity and mass-preserving external boundary conditions for 2D-doubly-parabolic model (3); • the introduction of mass-preserving and positivity-preserving permeability conditions at the interfaces between 2D and 1D parabolic models, see paragraph 3.2.2; • the introduction of mass-preserving and positivity-preserving permeability conditions at the interfaces between 2D-fully-parabolic model and 1Dhyperbolic-parabolic model, see paragraph 3.2.3. The plan of the paper is as follows. In section 2 we describe the biological framework that inspired our study, while in section 3 we introduce the mathematical formulation of biologically inspired models and we introduce the adopted model. Section 4 is devoted to the numerical techniques used to approximate the problem and in section 5 some numerical tests showing the qualitative behavior of cells in the designed environment are presented. Finally, in section 6 a discussion on the results and the future developments of our work is presented.

Biological framework
The control of immune cells migration and interaction with tumor cells living inside the chambers of the microfluidic chip, represent a new and attractive approach for the clinical management of tumor deseases. Furthermore, in the chip environment also drug testing can be exploited. Then, the quantitative assessment of immune cell migration ability to recognize and attack the tumor cells for each patient could provide a new potential parameter predictive of patient outcomes in the future. Migrating cells respond to complex chemical stimuli (as mixture of growth factors, cytokines and chemokines) representing a source of chemoattractants. These chemoattractants, through the interaction with their receptors allow cells to acquire a polarized morphology and to perform the action of immunosorveillance. The development of lab-on-chip technologies made it possible to realize a reproducible tailoring of the cellular microenvironment, thus allowing the continuous monitoring of experiments and accurate control of experimental parameters. Recently, the development of microengineering has given the possibility to realize culturing of multiple cell types and made it possible to observe cell-cell interactions and to transpose in vivo studies to a second generation of in vitro smart environments. The main advantages of this new technological tool are a close control over local experimental conditions and lower costs with respect to the use of animals in laboratory experiments for efficacy and toxicity testing. Some results obtained with on-chip experiments are presented in [1,7,25,29,38]. Regarding the structure of microfluidic devices, they are designed to allow chemical and physical contacts between tumor cells and non-adherent immune cells (i.e. murine splenocytes or human peripheral blood mononuclear cells). The microfluidic co-culture platforms are fabricated in polydimethylsiloxane (PDMS, Silgard 184), a biocompatible optically transparent silicone elastomer. In particular, here we refer to the experiment of two main culture chambers (a tumor and an immune cell compartment) connected via narrow capillary migration micro-channels having, respectively, width, length and height of 12µm, 500µm and 10µm. The cross-sectional dimensions of culture compartments are 1mm (width) × 100µm (height). For the development of the mathematical modelling explained in the next section we remark that we neglect the third dimension, thus we consider corridors and chambers as 2D objects. Two populations, immune and cancer cells, are introduced in two separate chambers; in particular, the immune cells are in the right chamber and the cancer cells in the left one (target chamber). The microchannels connect the two areas and allow the chemical diffusion and the migration. The culture medium is neutral, thus meaning that no exogenous substance are introduced. Mainly, the dynamics observed is the migration of immune cells from the right to the left in order to attack the tumor cells. The laboratory experiments are made in the context of immune competence vs. immunodeficiency, i.e. in a healthy or defective immune system. Indeed, the complex interactions, including cell-cell contacts, between cancer cells and immune system, which acts by limiting or suppressing tumor progression, is crucial in the tumor growth and invasion process.

Mathematical framework
Nowadays, mathematical analysis of biological phenomena has become an important tool to explore complex processes, and to detect mechanisms that might not be evident to the experimenters. Although a mathematical model cannot replace a real experiment, it may represent a support tool to explain acquired biological data and it may consent to gain a deeper understanding of the interactions between cancer cells and immune system. More generally, mathematical models can describe a  [38] edited by AAAS.  broad variety of biological phenomena, including cell dynamics and cancer [2], [14], [33], [26], [10]. The movement of bacteria under the effect of a chemical substance has been widely studied in the last decades, and numerous mathematical models have been proposed. As shown in [27], chemotaxis is decisive in biological processes. For instance, the formation of cells aggregations (amoebae, bacteria, etc) occurs during the response of the different species to the change of the chemical gradients in the environment. Moreover it is possible to describe this biological phenomenon at different scales. For example, by considering the population density as a whole, it is possible to obtain macroscopic models of partial differential equations.
In this paper, in order to describe the dynamics of cells in the 2D chambers we use a KS-like model, while in the microchannels we compare the behavior between two different modelization: 1D KS-like model and 1D GA-like model. The modelling here applied is described in the next subsection 3.1.
3.1. The model. Here we introduce a mathematical model that aims at describing the behavior of two populations of cells cohexisting together: tumoral cells T and immune cells(macrophage) M . We underline that the setting here considered can be make more complex with the introduction of a greater number of cell species and with the presence of an exogenous substance in the environment. The model consists of a reaction-diffusion system with chemotaxis, that it is able to describe birth and death processes, interaction with chemoattractant, interaction and competition between different cell species. The microfluidic chip is schematized as a network of channels connecting two boxes (the microfluidic chambers), then, following the ideas in [5], ad hoc transmission conditions were introduced to ensure the mass conservation. The parameters of the model, such as the velocity of different cell populations, the turning rates and the decay rates will be calibrated with observed data.
Cancer cells T produce chemical signal, called ϕ, activating the immune response of M and influencing their behavior. Moreover, we take into account the presence of cytokines ω (produced by M ), acting as a chemical killer of cancer cells. Then, the model describing the dynamics of the two cell species and the diffusion of the chemoattractant in the 2D chambers reads as: and we need to assign suitable initial conditions and boundary conditions for the cells and the chemoattractant concentrations that will be specified in the next paragraphs. In particular, the system above describes the following situation: tumor cells T produce a chemical substance ϕ attracting immune cells M and enabling them to recognize and interact with tumor cells. Immune cells also produce a chemical substance ω which makes the immune cells able to migrate towards the tumor cells. Therefore, in the first equation of the system (3), besides the diffusion term, we can find −λ T (ω)T representing the tumor suppression operated by immune cells.
In the second equation, in addition to the diffusion term we have the chemotactic term f = χ(M, ϕ)∇ϕ due to the presence of the chemical substance ϕ produced by the tumor. We remark that both in the first and in the second equation we include a term −k T (t)T and −k M (t)u taking into account the possibility of drug administration, with the functions k T and k M having an exponential decay in time: We also need to introduce the functions: representing, respectively, the chemotactic sensitivity of immune cells and the decay rate of cancer cells under the action of immune cells. Note that k ω is the killing efficiency of immune cells, k 1 represents the cellular drift velocity, while k 2 is the receptor dissociation constant, which says how many molecules are necessary to bind the receptors. We mainly refer to [27] for the values of the parameters k 1 , k 2 , γ, and all the parameters are reported in Table 1. Now, in order to describe the dynamics of cells in the microchannels connecting the two boxes, we introduce the following 1D models for the dynamics. To this aim, we consider two possible approaches to observe a different dynamics in the channels: • if we assign the 1D doubly-parabolic model on each channel, we have onedimensional version of system (3), where the superscript c indicates that we consider all the quantitites in the channels: • if we consider the 1D hyperbolic-parabolic model on each channel, we have the following system: We remark that, for the hyperbolic-parabolic system (7) we also need to assign initial and boundary conditions for the flux v.
For the sake of simplicity, we write the 2D model (3) as a general 2D-doublyparabolic system with source term as: with u the density of individuals and φ the density of chemoattractant. From now on, the two components of the drift term f will be indicated as: In the mono-dimensional channel we rewrite the 1D-doubly-parabolic system (6) in a more general form: and the 1D hyperbolic-parabolic system (7) with source term rewrites as: Table 1 are reported the parameters of the problem. The systems above have to be complemented with initial conditions for the unknowns u, v, φ, assumed to be smooth; initial data will be specified in paragraph 4.2.3. On the boundary we consider for all the quantities homogeneous Neumann conditions, so that we are assuming no-flux boundary conditions.
Monotonicity conditions. We also mention at this point that this model has an analytical monotonicity criteria. For linear convection term f c = au, and linear source term g = bu the criteria a λ − b ≤ 1, must be satisfied in order for the quantity u to be non-negative. Otherwise we would have negative u which would lead to unphysical solutions.
In regards to our model, that would mean we have for the immune cell density M the monotonicity condition and for tumor cell density T : to be verified in the computational domain in order to ensure non-negative solutions.

Remark 1.
We remark that the no-flux conditions boundary conditions used in our simulations are needed to have the mass-conservation of all the quantities. However, they are not realistic, since in the laboratory experiment there is an inflow of cells from the outer boundaries. In our future developments we will extend the no-flux boundary conditions to more general ones.
For our implementation we want to model this PDE in two domains. A simplified schematization of the bounded surface where experiment is performed is reported in Fig. 2. We have two microfluidic chambers of the same size, one on the left and the other on the right, defined, respectively, as Ω l = [0, L x ] × [0, L y ] and Ω r := [L x + L, 2L x + L] × [0, L y ] they are connected by microchannels, each of them schematized for simplicity as a line I = [0, L]. Thus, the link between the box on the left and the corridor is schematized as a junction (node 1L) and analogously the link between the corridor and the box on the right, as node 2L. The two junctions are not really a single point, thus they are parametrized as an interval for node 1L and node 2L, namely [a 1 , b 1 ] of length σ := b 1 − a 1 . We remark that for the sake of simplicity, the numerical treatment is developed for a simplest geometry composed by 2D chambers connected through a single 1D channel. The extension to multiple 1D channels is done in paragraph 4.2.2.

3.2.
Outer and inner boundary conditions for the models with source term g = 0. From now on, in order to use the mass conservation argument at the outer boundaries of the 2D domain and at the inner interface between 2D and 1D domains, we study a simplified version of the models (8)- (10) and (8)-(11) putting the source term g equal to zero: In particular, we have to prescribe the flux conservation at the inner boundaries for the 2D-1D parabolic case (14) and for the 2D parabolic-1D hyperbolic case (15), since we cannot loose nor gain any cells during the passage through a node. While keeping most of the boundary conditions, we must change them at the interface node. In the 2D left box Ω l , the position of node 1L is at x = L y , y ∈ [a 1 , b 1 ] and for the 1D domain represented by the channel node 1L is placed at x = 0, see Fig. 2.
We shall note at this point that although we work in the following with the general systems (8), (10) and (11), the same results can be applied to the complete models (3),(6) and (7). 3.2.1. Boundary conditions for the 2D doubly-parabolic model (14)-a). Considering that our model describes the migration of cells by both diffusion and chemoattractant effects, physically speaking the mass of cells and the chemoattractant must be preserved in absence of creation and destruction of cells. For φ in (14)-a), by using the divergence theorem, we can write: We assume no-flux condition for the chemoattractant in order to preserve its mass in absence of source terms. Since we defined the source term f as a product function of ▽φ, we get the equivalent condition f (t, x, y) n| δΩ = 0. Note that for u the same condition holds: and these boundary conditions guarantee mass-conservation.
The same approach gives no-flux condition for both u and φ in the right chamber Ω r (and for the complete model for T and ω as well).

3.2.2.
Interface between 2D-1D models in (14). Here we prescribe the conservation of the mass between the left box and the corridor (node 1L in Fig. 2). The conservation condition reads as: and it rewrites as: by using the divergence theorem in the first integral. With our analytical boundary conditions (17), the integral vanishes except at the boundary where the node is positioned.
We remark that attention has to be paid with n being the outer normal of the domain. We have: and, thanks to the boundary conditions (16) and (17) some terms cancel in equation above, thus we get the condition: Now we impose Kedem-Katchalsky (KK) [22] conditions describing the conservation of the flux through a node (see also [34] for numerical treatment of these conditions). In particular, at the interface between left chamber and channels we have (on the left of node 1L in Fig. 2): and on the right of node 1L we have: Thanks to conditions (19) and (20) we are guaranteed to have the flux conservation (18); and we will use such conditions to obtain numerical boundary conditions for the boundary values at the nodes on both sides, as shown in Section 4 in paragraph 4.1.2.

3.2.3.
Interface between 2D-1D models in (15). In this section we describe the combination of 2D parabolic-1D hyperbolic model in order to describe the dynamics with a hyperbolic model (11) in the one-dimensional domain represented by microchannels. Further care has to be made in order to keep some important properties which ensure consistency and non-negativity of numerical solutions when connecting both models. Now the transmission condition for the switch from Ω l to I = [0, L] are derived in this case. For the mass conservation we impose the condition: Note that in the above formula we have v (L, t) = 0 because we are looking at left interface (node 1L). Then, we finally get: Now we impose the KK-condition at the interface: and then (21) reads as:

Numerical approximation
Here we describe the numerical approximation of the adopted models, 2D-doublyparabolic, 1D-doubly-parabolic and 1D-hyperbolic-parabolic. We define equispaced x i := i△x, t n := n△t and y j := j△y with △x, △y, △t > 0 and i = 0, . . . , N x + 1, j = 0, . . . , N y + 1; for the channel [0, L] we discretize it as x i = i△x, with i = 0, . . . , N . For a more structured presentation, we introduce the operators We use a first order explicit finite difference method in time and a second order central method for the approximation of the diffusion term in space. For the chemoattractant term we use a finite difference scheme in space as will be specified in the sequel. We remark that using a purely explicit methods implies restrictions on the mesh grid spacing and time step size to ensure stability due to the Péclet number criterion. However, to prevent strong restrictions on the mesh grid in the case of dominant advection regime, we introduce artificial viscosity, which leads to less restrictive condition when diffusion is small, but it decreases the order of the scheme. In the sequel we always assume to have the mesh grid small enough, thus neglecting the artificial viscosity term, which will be addressed in the numerical approximation of the model in section 4.
Remark 2. Note that special attention has to be paid also to the source term g (x, y, t, u). Although in a simple explicit method one can evaluate the function at each time step n at mesh grid point (i, j), the function itself can induce stiffness, enforcing small time steps. To overcome this issue, implicit methods can be used such as the Crank-Nicolson-method. However, here we work with pure explicit methods to make an easier presentation of the schemes and we will address the implicit method for these models in the sequel, see paragraph 4.2.3.
Another issue is the choice of the right boundary conditions which should reflect the qualitative attributes of the analytical model. In absence of source terms, the mass of cells and chemical substances are preserved. In order to make a numerical verification of this property, we considered the numerical approximation at the interface between 1D-1D models. In more detail, choosing standard boundary conditions by simply discretizing Neumann boundary conditions with a finite difference scheme, the mass will not be preserved over time, see Fig. 3. In particular, in Fig. 3 a comparison between mass-preserving and usual finite difference boundary condition is performed, for the 1D-doubly parabolic case on both sides of the interface (on the left) and for the 1D-doubly-parabolic-1D-hyperbolic-parabolic interface. From this 1D numerical example it is evident the necessity to develop modified boundary conditions which are consistent and preserve the mass correctly.
Mass-preserving and positivity-preserving numerical approximation will be developed in the present section. In the following we will neglect the label c to make the reading easier and make distinction only when necessary.  . On the left: evolution of total mass for 1D-doubly-parabolic model with standard(finite difference) and mass-preserving boundary conditions. On the right: evolution of total mass for 1D-hyperbolic-parabolic model with standard (finite difference) and mass-preserving boundary conditions.

4.1.
The parabolic-parabolic case. Here we propose a numerical scheme for the doubly-parabolic systems (8) and (10).
For the discretization of equations in 2D system (8) in the interior points of the domain, i.e. for i = 1, . . . , N x , j = 1, . . . , N y , we define an explicit in time finite difference discretization both for u and φ: consistent approximation for div (f ), i.e. a second order central in space finite difference: and analogously for △ n j (f y i ). Note that the function f x = χ(u)φ x can be discretized with f x,n i = χ(u n i,j )(φ n x ) i,j with (φ n x ) i,j an appropriate second order approximation of φ x : For the chemoattractant we have the approximation scheme: +au n i,j − bφ n i,j . For 1D system (10) in the interior points of the channel we apply the same explicit in time finite difference scheme both for u and φ: with △ n i (f ), as above, a second order central in space finite difference: where the 1D version of (24) for the approximation of φ x in f is used.
CFL condition. By using the Von-Neumann stability analysis we obtain the following stability criterias (CFL-condition).
These restrictions for the step size and the mesh grid size can be avoided by using implicit methods, but this would increase computational cost because of the necessity of solving a non-linear equation system at each iteration.
Remark 3. Since we are dealing with an explicit method, we can calculate the values for the next time step n + 1 by using solely the values of the previous time step n; we remark that an implicit approximation can be applied, with the use of Crank-Nicolson (CN) method, in order to increase the accuracy to second order and avoid restriction of △t and △x due to the CFL-Condition which do not arise with the CN-method. Indeed, we remark that the implemented algorithm for simulations described in paragraph 4.2.3 is based on the CN-method in time.
For the two-dimensional system here we report the numerical scheme. The numerical method in the interior points of the 2D domain for the cell density reads as: The numerical method in the interior points of the 2D domain for the chemoattractant reads as: In the following we present the discretization of the boundary and transmission conditions to complete the numerical schemes. 4.1.1. Discretization of the boundary conditions for the doubly-parabolic problem. Now, in order to complete our numerical scheme, we need to discretize the boundary conditions to obtain values for the boundary on each domain for the time step n + 1. Since a qualitative characteristic of this model is the preservation of total mass, we want our numerical model to preserve mass at each time step. To this aim, we have to choose discrete boundary conditions that both are consistent with the analytical boundary conditions and preserve the mass in the numerical method. We remark that we present the computations without source term g and we will add it in the sequel to complete the equations.
Boundary conditions for the density of individuals u. The mass conservation over time on Ω l reads as: Now, applying a quadrature rule for the numerical integration: we need to ensure that For the numerical integration different quadrature formulas can be used. Since we want to use constant space-steps and want to obtain mass-preserving boundary conditions for the numerical methods, closed Newton-Cotes methods are suitable. In particular, we use the trapezoidal rule which introduces an integration error of O △x 2 . For the one-dimensional trapezoidal rule with a function z : and for the two-dimensional trapezoidal rule with function z : R 2 −→ R we have: Imposing the equality I n+1 − I n = 0 in the 1D case gives: Using the numerical scheme (26) for u n+1 i for i = 1, . . . , N we get: Then we obtain: Then, applying (28) we obtain in the above formula: We can now compute the values for both u n+1 0 and u n+1 N +1 so that the term equals to zero. By collecting values from nearby stencils together (otherwise we obtain an error of O(△x) which can be verified by Taylor expansion), we obtain the following conditions at the outer boundaries of 1D domain: and In the interior points of 1D domain we have the numerical scheme: Then we can state the following result.
Proposition 1. The scheme (41) endowed with boundary conditions (39) and (40) is mass-preserving by construction, since it is obtained imposing I n+1 − I n = 0, as shown above. Moreover, the scheme, obtained with the integral method above, is second order in space up to the boundaries since it can be equivalently obtained using the following second-order approximation of the first derivative including a ghost cell: Finally, the scheme is also positivity-preserving under the parabolic CFL condition.
Now we compute u n+1 0 directly by using the second-order centered numerical scheme and replace the ghost value u −1 from the discretization of condition (42). While this works well when f = 0, the same does not happen for f = 0, thus making the approach with the discrete integral equation still necessary. Futhermore, by using a different numerical integration scheme, we can achieve different masspreserving boundary conditions of higher order. Using the mass-preserving property argument, we compute boundary conditions for the corners and top and bottom boundaries of the 2D domain Ω l for f = 0. By applying them with the numerical method (31) into I n+1 − I n = 0, we get the expression: since the terms in u cancel. By choosing again the central in space second order finite difference scheme (23) for div (f ), we get Now we can distribute the remaining values to the boundary values in the same way we did for the 1D-parabolic case. Therefore, we obtain the following masspreserving boundary conditions: − △t △y f y,n Nx+1,Ny+1 + f y,n Nx+1,Ny and for the top and bottom boundaries we have: Boundary conditions for the density of chemoattractant φ.
For the computation of the conditions at the outer boudaries for the chemoattractant φ c in the 1D-doubly parabolic model we proceed as above, but neglecting the source term a c u − b c φ c to obtain boundary conditions that are mass-preserving. By doing so, we achieve the following second-order accurate and mass and positivity preserving boundary conditions for the chemoattractant: The parabolic equation in the interior points is solved using a finite differences scheme in space and an explicit method in time: In a similar way we can extend the numerical boundary conditions for the 2Dparabolic model.

Discretization of the transmission conditions for the doubly-parabolic case.
Since for K = 0 we would achieve the same analytical boundary conditions for separate domains (17) and (16), we follow the same approach by using ghost values as in (42) in the numerical scheme to obtain the boundary conditions, since in such a way mass preserving and positivity preserving boundary condition are achieved. By using the approximation formula (28) in the condition (19) on the left of node 1L we have: Then we have: Nx+1,j and we get: for j = j a1 , . . . , j b1 .
Moreover, using (28) in (20), we can write: y, t)) dy + f n 0 and then we obtain: and we finally get the formula: We now use the Ansatz to apply the ghost values into the numerical scheme without specific chemotactic approximation (41) and (31), and use the discrete integral equation to determine the chemotactic term discretization. Because we now not only need to conserve the mass in each domain, but in both connected ones, the must use the expanded discrete integral equation to compute the total mass over both domains.
Plugging the ghost values (50) and (52), respectively, into the numerical schemes (31) and (41), we get the conditions at the interface (node 1L): In particular, the conservation of the discrete total mass reads as: and now, applying the conditions (53) and (54) with the other boundary conditions (44) and (39) we get: and obtain the following transmission conditions (56) for j = j a1 , . . . , j b1 . The integral expression of the density u n+1 0 in (56) can be expressed with a numerical quadrature form, such as the trapezoidal rule, as in (36). Proceeding analogously as above, this approach leads to mass-preserving and positivitypreserving transmission conditions for the chemoattractant φ as well. In particular, we have at the first and last endpoint, respectively: We have finally developed a complete numerical scheme to treat doubly-parabolic partial differential equations systems in two domains, 1D and 2D, connected through a node, which ensures the mass conservation and the positivity as the original PDE.

4.2.
The hyperbolic-parabolic case. The second order AHO scheme on a line was introduced in [28] for the 1D hyperbolic system (2). Here, considering the presence of the source term g on the right hand side of the equation for the density of cells, the AHO scheme reads as: with mass-preserving boundary conditions (including the additional source term g) at the external boundaries. We remark that for the hyperbolic-parabolic model not only mass must be preserved as the in the fully-parabolic model, but also the flux v needs to converge towards the steady state v = 0. Since here we have the 1D domain connected at both the endpoints we do not need to use numerical boundary conditions for the outer boundaries. However, for the details and the description of the AHO numerical scheme at the outer boundaries, see [28]. For this reason we use the so called AHO (Asymptotic Higher Order) schemes (see [6] for the study of AHO scheme at interfaces including mass-preserving transmission conditions) with source term g for which the approximation of the stationary solutions is up to third order of accuracy and it converges towards a numerical solution with v = 0, while preserving the mass.

Discretization of transmission conditions for the 2D-doubly-parabolic and
1D-hyperbolic-parabolic case. The first equation is the same as for the 2D-doublyparabolic and 1D-doubly-parabolic case. Hence we derive the same transmission condition for u n+1 Nx+1,j for j = j a1 , . . . , j b1 . For the flux, the transmission condition (22) gives us We remark that here we have two problems which do not occur in the previous model. In particular: we do not have a transmission condition formula for u n+1 0 and the formula for v n+1 0 is implicit, since it depends on values at time step t n+1 . In order to solve these issues, we will use once again the equivalence between discrete total masses to obtain mass-preserving computation formula for u n+1 0 and, since v n+1 0 only depends on u n+1 0 and u n+1 Nx+1,j , which can be obtained explicity from computed values at time step t n , we can calculate v n+1 0 explicity as well. Then, imposing that: and we finally obtain the transmission condition with source term g : Proposition 2. The complete numerical scheme derived for the 2D-doubly-parabolic-1D-hyperbolic-parabolic model has the feature to be mass-preserving across the transmission conditions. Note that for the chemoattractant equation is the same as for the 1D-doubly-parabolic and 2D-doubly-parabolic case. Hence, the numerical schemes (47) and (32) with boundary conditions (45) and (59) can be used.

Multiple channels.
In the previous paragraphs we have connected the twodimensional domain Ω l with a single one-dimensional channel I at (L x , y) ∈ Ω l with y ∈ [a 1 , b 1 ], and j a1 and j b1 , the positions of the endpoints of the corridor on the numerical grid. Of course this can be extended to more channels.

Implemented algorithm.
Before presenting the numerical tests in the next section 5, we detail the approximation scheme for the density u, also including the source term g, implemented to solve the problem in the 2D-1D domain. As underlined before, it is necessary to use implicit schemes to consider the presence of stiff source terms. For this reason, for the approximation of the time derivatives we use the Crank-Nicolson method on the diffusion and source term, which is a second order implicit method and the explicit central method for the convection term. Because of the explicit term, we have numerical restrictions on the mesh grid and time step. Furthermore, as discussed previously, we introduce artificial viscosity to avoid oscillations due to not suitable mesh grid size in dominant convection regime, which is often the case in chemotaxis models. The implicit-explicit numerical method used to compute the solutions for the density u in (8) inside the 2D domain Ω l is: artificial viscosity , with θ n i,j := χ(u n i,j , ϕ n i,k )|∇ϕ n i,j | for i = 1, . . . , N x , j = 1, . . . , N y . As can be seen, the function θ used for the artificial viscosity is almost identical to f with the exception of using the absolute value of ∇ϕ. By using this, we increase artificial viscosity only where the gradient of the chemoattractant increases. This reduces the restriction on the meshgrid due to the condition induced by the cell Péclet number. The numerical transmission condition on the left of node 1L (i = N x + 1, j = j a1 , . . . , j b1 ) is: Nx+1,j additional term for transmission condition .
The role of KK coefficient K in the positivity of (64) is discussed in Remark 5.
For the corners we use the following boundary conditions:  .
Similarly, for the chemoattractant φ we have the implicit-explicit scheme in the interior points of the 2D domain: and for the boundaries and the corners the numerical schemes for φ are, respectively, (48) and (49).

Remark 4.
If we consider two-dimensional domain Ω r connected to the right endpoint of the one-dimensional corridor I, the complete numerical scheme for the left domain Ω l described above can be considered.
The main difference is that the transmission conditions at the interface between the box and the channel (the left for the box Ω l and the right for the corridor) are reversed to the left for the corridor and the right for the box Ω r . In the numerical scheme, the only change affects the channel I, where we have transmission conditions also for u n N +1 (resp. v n N +1 ). The same boundary condition can be used without transmission conditions, with only the additional term derived from the KK-condition and it must be added as well for u n 0 (resp. v n 0 ).
For the computation of solutions on the one-dimensional channel I, we have two different approximations depending on the choice of the model we assign on it. If we solve the doubly-parabolic problem (10), the approximation scheme used is the Crank-Nicolson scheme, as above: with the transmission condition on the left of node 1L (i = 0) given by: same as for BC without transmission condition If, instead, we need to solve the hyperbolic-parabolic problem (11), an implicit version of the second order AHO scheme (60) is used. Indeed, because of the different scales of the parameters, instabilities in the corridors can occur if the stability condition of the AHO scheme is not satisfied. In particular, inside the channels we use the scheme: , endowed with the following implicit version of transmission condition (61) : and analogously for the condition (62).
Remark 5. Note that, in order to ensure the positivity of the quantities in the above formulas deriving from the KK conditions, i.e. (64) for the 2D domain and (69) or (70) for the 1D domain, we also need to take care of the ratio between the KK coefficient K and the space discretization steps. In particular, for (64) and (70) one needs to ensure that K △t △x and, respectively, K △t △x △y is not too big in order to damp possible high obscillations produced by the term in parenthesis. Similarly, in (69) we need to check that K △t △x σ is small in order to prevent the growing of negative term.
Moreover, as previously discussed, we need to check that the numerical monotonicity conditions are satisfied: in the computational domain in order to ensure non-negative solutions.
For the sake of completeness, we underline that at each time step a non-linear equation system must be solved, for which Newton-Krylov-subspace methods [24] can be used which take advantage of the mostly sparse structure of the jacobian matrix.

Numerical tests and results
This section is devoted to the presentation of the numerical tests and the parameters of the problem are reported in Table 1. Our aim is to show the ability of the simulation algorithm based on the model (3)-(7) to reproduce the qualitative behavior of the two population sharing the same habitat as observed in the videos of laboratory experiments. We remark that we decided to do numerical simulations of the chip geometry assigning the 1D-hyperbolic-parabolic model on channels since it seems more realistic. However, a numerical test on the behavior of the model (3)-(6), with the doublyparabolic model on channels, is provided in the last Example 4.
Example 1. Before we numerically simulate the laboratory experiment with the algorithm, we conduct a simple numerical test in order to prove its accuracy. We assumed the following setting: a left squared chamber Ω l with one corridor positioned in the middle and only one cell family with initial distribution u(x, y, 0) = 5e − 1 2 ((x−0.5) 2 +(y−0.5) 2 ) . Since we do not have any analytical solution for this problem, we choose dt and dx = dy small enough to obtain reasonable error estimations. In this case we use dt = 10 −4 and dx = dy = 5 × 10 −4 for the approximation u e at time t = 100 and calculate the error as the quantity u e − u approx in L 1 -norm. In order to confirm the order of our scheme, we use a log-log-plot with constant and small enough dt (resp. dx), and decreasing dx (resp. dt). As shown in Figures 4 and 5 the time order and space order equals to line with slope 2 in the log-log plot which corresponds to our scheme of order 2 in space and time.     (x, y) ∈ Ω l is chosen as: (73) T (x, y, 0) = 5e − 1 2 (x 2 +y 2 ) + 5e − 1 2 (x 2 +(y−5) 2 ) + 5e − 1 2 (x 2 +(y−10) 2 ) , whereas in the corridors and the right chamber no tumor cells are present. For the immune cells distribution on the chip for (x, y) ∈ Ω r we assign: whereas no immune cells are present in the left chamber nor in the corridors. For the chemoattractants we set a constant initial density: ω(x, y, 0) = 0 (all domains) and ϕ(x, y, 0) = 2 for x, y ∈ Ω l \ {L x × [a 2 , b 2 ]}, ϕ(x, y, 0) = 0 for x, y ∈ Ω r − [a 2 , b 2 ] and a linear decreasing in space initial value ϕ(x, y, 0) = −0.01x + 5 in correspondence of nodes 2L − 2R, namely at (L x , y) and (L x + L, y) for y ∈ [a 2 , b 2 ].
For this simulation test we choose the parameters for each domain as given in Table  1 for both the chambers and for all the corridors we used the same parameters. The numerical method implemented is listed in paragraph 4.2.3; for the 1D channels the AHO-Scheme (70) is implemented since we are considering the hyperbolicparabolic model. The discretization grid has time step size △t = 10 −3 and space size △x = △y = 0.25.    In Figures 6, 7, 8 and 9 we can see the density of the tumor cells T and immune cells u for different times t = 0, t = 5, t = 10 and t = 100 accordingly. Note that at time t = 0 tumor cells are present in the left chamber only and immune cells are present in the right chamber only. Since no chemoattractant is present at the initial time t = 0, both cells diffuse around their chamber and slowly entering the corridors while creating chemoattractant ϕ and ω. But already at time t = 5 we notice that in the middle of the left chamber the tumor cells are getting pushed towards the buttom and top of their chambers. Indeed, tumor cell densities T slowly diffuse around the left chamber but accumulate around the bottom and top and partly enter the corridors. This is due to the fact that the chemoattractant ϕ produced by cancer cells (the annexin) induce a migration of the immune cells M towards the tumor cells T causing a higher migration towards the center of the left chamber where the initial distribution of tumor cells was closest to the chambers. For t = 10 in Fig. 8 and t = 100 in Fig. 9 we see that, since most tumor cells are accumulating on the top, they manage to diffuse through the nearest corridor on the top into the right chamber; on the other hand, the immune cells continue to migrate towards the highest concentration of chemoattractant ϕ, which is where the tumor cells T concentrate. Especially at time t = 100, we can see that the quantity of tumor cells has dramatically decreased compared to the immune cells caused by the action of chemokine ω.        We can clearly see in Fig. 11 compared to Fig. 7, that the immune cells u are much more massively moving towards the left chamber due to the chemoattractant ϕ, which causes a slightly higher concentration in the left chamber. But due to the diffusion of the chemoattractant ϕ and the creation of more chemoattractant from the tumor cells, the graphs of both Examples 2 and 3 are getting similar during the time evolution until a difference is no more noticeable, around time t = 100.
Example 4. In this last Example, we tested the 1D-doubly-parabolic model on channels and compared it with the hyperbolic-parabolic model used in the previous Examples. In Figures 14-16 we assigned the 2D-doubly-parabolic model which uses the parabolic partial differential equation to describe the movement in the 1D channels. By using the same initial data as for the other Examples, we notice that for time t = 100, the doubly-parabolic model in Fig. 16 seems to have a similar pattern as for the hyperbolic-parabolic model depicted in Fig. 9 and the hyperbolic-parabolic model with stronger chemotaxis in Fig. 13, but the scale differs a lot between these models. Whereas we have for the tumor cells T a maximum concentration of 10 −5 for the hyperbolic-parabolic model, and 10 −7 for the hyperbolic-parabolic mode with stronger chemotaxis, we see clearly that for the doubly-parabolic model, the concentration of the tumor cells T is of the order of 10 −2 . This is due to the much slower movement of the immune cells through the corridors. This also explains the much higher concentration of the chemoattractant φ because of the much higher concentration of T compared to the other models.
In the following Figure 17, we represent the density of tumor cells and immune cells depicted in Figures 6-9 as individuals, by randomly placing them according to their density. The higher the density at a given point, the more cells will be distributed randomly around that area. If the density is lower than a chosen threshold in a certain point, no cells will be represented around it.  (c) Visualization for time t=50. Figure 17. Visualization of immune cells (blue dots) and tumor cells(red squares) for time t=0, t=5 and t=50 by using the density of each quantity and representing them as cells.

Conclusion and future perspectives
The principal feature of the present work has been the development of a simulation tool to describe cell movements and interactions inside microfluidic chip environment. Our study focused on both the modelling and the numerical point of view. Indeed, schematizing the chip geometry as two 2D-boxes connected by a network of 1D-channels, the main issues were: • the introduction of mass-preserving conditions involving the balancing of incoming and outgoing fluxes passing through interfaces between 2D and 1D domains; • the development of mass-preserving numerical schemes at the boundaries of 2D domain and mass-preserving transmission conditions at the 2D-1D interfaces.
Furthermore, from the modelling point of view, we studied the dynamics in the channels in case of doubly-parabolic model and hyperbolic-parabolic model. Since we obtained comparable asymptotic states, we decided to apply the hyperbolicparabolic model in order to have finite speed of propagation in the channels which seems to be more realistic. In this framework, having in mind the laboratory experiments on chip described in section 2, it was possible to simulate the chip environment with two species of living cell moving in it. Moreover, we remark that we can simulate more complicated situations where more than two cell species are present. As a further development of the present study, we will work on the calibration of the model against experimental data.