Mass-Preserving Approximation of a Chemotaxis Multi-Domain Transmission Model for Microfluidic Chips

Abstract

The present work is inspired by the recent developments in laboratory experiments made on chips, where the 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 is 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 a 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.

1. Introduction

The aim of the present work is to study both the modelling and numerical approximation of a chemotaxis-reaction-diffusion mathematical system describing the qualitative behavior of different cell species living in a confined environment. This work is inspired by laboratory experiments made on microfluidic chip [1], where some populations cohexist and interact. In recent years, there has been the development of a new approach to biological studies aimed at reconstructing organs and complex biological processes on a chip [2]. 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 the laboratory setting of the experiment in microfluidic chips [1,2,3]—see also the short description of the experiments reported in Section 2.1.1—we introduce a model mimicking the interactions between two cell populations—namely, immune and cancer cells.

The mathematical model, proposed in Section 2.2, is a reaction-diffusion system with chemotaxis and describes birth and death processes, the migration of immune cells driven by chemical signals produced by tumor cells, and interactions between different cell species. We underline that, since the chemical gradients are not measured experimentally, by using the simulation algorithm based on the mathematical model proposed here, the chemical concentration gradients in the chip can be obtained by solving the inverse problem of minimizing the residuals between the measured trajectories and the simulated ones; see also the discussion in Section 5 about the future developments of our work. From a mathematical point of view, we follow the framework of the classical macroscopic models of chemotaxis; see, for instance, [4], 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 the 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 [4] 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 (model and numerical algorithm) could be easily extended to the third dimension and we remark that, analogously to the 2D-1D case here considered, mass-preserving and positivity preserving numerical conditions at the external boundaries and at the interfaces between 3D and channels still hold.

We consider the microchannels connecting the 2D chambers as 1D lines for modelling and computational reasons, as explained in Section 1.1. In order to model the dynamics on the microchannels, we choose two different approaches: we can assign a 1D version of the doubly parabolic model used in the chambers; otherwise, we can assign a model derived from a 1D-GA model [5], 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 [6,7,8,9,10].

1.1. The Geometry of the Microfluidic Chip and of the Related Computational Domain

The microfluidic chip is represented as a network of channels connecting two boxes (the microfluidic chambers); see Figure 1 and a schematic picture of the related computational domain is depicted in Figure 2. 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 with, respectively, width and length of $12\mathsf{\mu }$m and $500\mathsf{\mu }$m. Moreover, the channels height is of $10\mathsf{\mu }$m; however, since in the video footage the experiments is recorded at a fixed height, the third spatial dimension in our framework is neglected. The cross-sectional dimensions of culture chambers are 1 mm (width) ×$100\mathsf{\mu }$m (height).

A simplified schematization of the bounded surface where the experiment is performed is reported in Figure 2. We have two microfluidic chambers of the same size, one on the left and the other on the right, defined, respectively, as ${\Omega }_l=[0,L_x]\times [0,L_y]$ and ${\Omega }_r:=[L_x+L,2L_x+L]\times [0,L_y]$. They are connected by microchannels, each of them schematized as rectangles $R=[0,L]\times [a,b]$. In order to ease the reading, we point out that in the sequel we approximate the rectangular microchannels as 1D intervals $C=[0,L]$ with zero thickness for the following reasons:

- 

modelling reason the width of microchannels (12 $\mathsf{\mu }$m) is comparable to the size of cells (for instance, immune cells measure about 8–10 $\mathsf{\mu }$m of diameter);

- 

computational reason to reduce the running time of the simulation algorithm, since otherwise we should consider a 2D meshgrid for each microchannel.

Then, 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, therefore they are parametrized as an interval for node $1L$ and node $2L$—namely, $[a,b]$ of length $\sigma :=b-a$.

We remark that for the sake of simplicity, the numerical treatment is developed for a simpler geometry composed by 2D chambers connected through a single 1D channel. The extension to multiple 1D channels is done in Section 3.2.2.

1.2. Original Contribution of the Present Paper

From the mathematical and numerical viewpoint, here we deal with a challenging issue arising in the 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 [4] of chemotaxis reads as:

$$\left\{\begin{array}{c}{u}_t=div\left(\nu \nabla u-\chi \left(\varphi \right)u\nabla \varphi \right)\hfill \\ {\varphi }_t=D\Delta \varphi +au-b\varphi ,\hfill \end{array}\right.$$

(1)

with u the density of individuals in the considered medium, $\nu$ the diffusion rate of the organism according to Fick’s Law, and $\varphi$ the density of the 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 $\chi$ 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 Section 2.2.1.

In the 1D microfluidic channels, we use the one-dimensional version of the KS-like model used in the chambers, but we also study the behavior of individuals when a hyperbolic-parabolic model, characterized by finite speed of propagation, is assigned. Such a hyperbolic-parabolic model, described in Section 2.2.1, is inspired by the Greeberg–Alt (GA) model, arising as a simple model for chemotaxis on a line:

$$\left\{\begin{array}{c}{\partial }_{t}u+{\partial }_{x}v=0,\hfill \\ {\partial }_{t}v+{\lambda }^2{\partial }_{x}u=-v+\chi \left(\varphi \right)u{\partial }_x\varphi ,\hfill \\ {\partial }_t\varphi =D{\partial }_{xx}\varphi +au-b\varphi .\hfill \end{array}\right.$$

(2)

Note that here v is the averaged flux. Let us underline that the flux v in model (2) corresponds to $v=-{\lambda }^2\nabla u+\chi \left(\varphi \right)u\nabla \varphi$ for the KS system.

We remark that here we do not consider the GA model on the 2D domain, since in this case the derivation of the monotonicity condition, easily computed in 1D, it is not straightforward, due to the oscillations brought by 2D wave equation. One possibility to overcome this problem should be to consider the macroscopic hyperbolic model proposed by Preziosi in [9,11]. Alternatively, a kinetic 2D model of chemotaxis and its numerical approximation is studied in [12].

This system was analytically studied on the whole line and on bounded intervals in [13], while an effective numerical approximation, the Asymptotic High Order (AHO) scheme, was introduced in [14]—see also [15,16]—and extended on networks with general boundary conditions in [17,18].

Here, in the numerical treatment for the computation of numerical solutions one has to take care of what happens at the interface when switching from 2D-doubly parabolic models to 1D-doubly parabolic or 1D-hyperbolic-parabolic ones and vice versa.

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 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 3. The approximation of doubly parabolic chemotaxis models for the 1D-KS model (1) on networks was already considered in [19]. 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 required to impose the continuity of the density of both cells u and chemoattractant $\varphi$, 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 [14] for a single line. In particular, the numerical treatment of the hyperbolic part of the system is based on the finite difference AHO scheme with the development of mass-preserving numerical scheme at outer boundaries, while the parabolic part is approximated by finite difference and Crank–Nicolson scheme.

In [17,18] the GA system was solved on networks, thus making it necessary to develop mass-preserving transmission conditions at inner nodes (interfaces) and suitable boundary conditions at outer nodes. However, rgw transmission conditions considered there involved the mass exchange only 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 the case of oscillating functions.

The numerical approximation of permeability Kedem–Katchalsky [20] conditions describing the conservation of the flux through a node was already considered in [21], 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 [22] for finite difference schemes and in [23] for discontinuous Galerkin methods. The numerical treatment of permeability conditions for chemotaxis problems was presented for the first time in [24] 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 the interfaces.

Therefore, to the best of our knowledge the present paper is the first numerical work where this new technique of switching the size of the domains and type of equations (parabolic vs. hyperbolic approach) is introduced, in order to develop mass-preserving and positivity preserving schemes. We point out that in the present paper the approximation method—based on finite difference schemes—involves the derivation of suitable zero-flux boundary conditions. The motivation for our framework based on finite difference methods seems to be a natural choice when dealing with hyperbolic models of transmission.

1.3. Main Contents and Plan of the Paper

In the present paper, a positivity-preserving and mass-preserving numerical discretization of Neumann boundary conditions at the corners and at the bottom and top boundaries of the 2D domain for a 2D-doubly parabolic reaction-diffusion problem are presented. Moreover, a positivity-preserving and mass-preserving numerical scheme at the interfaces of the network connecting the 2D chambers with the 1D channels (where the 1D-doubly 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 categorized 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 Section 2.3.2;
• the introduction of mass-preserving and positivity-preserving permeability conditions at the interfaces between the 2D-fully parabolic model and 1D-hyperbolic-parabolic model—see Section 2.3.3.

The plan of the paper is as follows. In Section 2 we describe the biological framework that inspired our study and we introduce the mathematical formulation of biologically inspired models and we introduce the adopted model. Section 3 is devoted to the numerical techniques used to approximate the problem and in Section 4 some numerical tests showing the qualitative behavior of cells in the designed environment are presented. Finally, in Section 5 a discussion on the results and the future developments of our work is presented.

2. Materials and Methods

The present Section is devoted to:

• the description of the biological framework and the laboratory experiment that inspired our work—see Section 2.1;
• the mathematical methods—see Section 2.2—bringing us to the development of a simulation algorithm designed for qualitatively reproducing the experimental observations.

2.1. 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 diseases. 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 a 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 immunosurveillance.

The development of lab-on-chip technologies has made it possible to realize a reproducible tailoring of the cellular microenvironment, thus allowing the continuous monitoring of experiments and the accurate control of experimental parameters. Recently, the development of microengineering has provided 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,2,25,26,27].

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.

2.1.1. Setting of the Laboratory Experiments

Here, we shortly describe the laboratory experiments inspiring our work (see [1]) and carried out in the microfluidic chip environment; see Figure 1. Two populations, immune cells of wild type and cancer cells, are initially put into two separate chambers and can have physical and chemical contact through the microchannels. In more detail, rgw cells are loaded into the reservoirs as follows: the left chambers are filled with about $2\times {10}^6$ human peripheral blood mononuclear cells and the right chambers with about $5\times {10}^4$ breast cancer cells, pre-treated or not with doxorubicin hydrochloride, all immersed in a suitable culture medium. Time-lapse recordings are performed over a period of 72 h (1 microphotograph every 2 min) by means of a microscope placed directly inside the CO${}_2$ incubator for the duration of the recording.

We consider two different scenarios:

• first scenario (treated case): before enabling cells to migrate, tumor cells are previously treated with a chemoterapy drug. Afterwards, we observe immune cells migrating towards the left chamber where the tumor remains confined, but expresses the chemical stimuli attracting immune cells. Mainly, the dynamics observed in this case is the migration of immune cells from the right to the left in order to attack the tumor cells.
• second scenario (untreated case): tumor cells migrate in the right chamber and proliferate. In this case, the tumor cells do not produce chemoattractant, thus immune cells move in the environment without recognizing and attacking tumor cells.

The culture medium in both cases is neutral, thus meaning that no exogenous substances are introduced. Our aim is then to build a simulation algorithm based on the mathematical model, which is able to reproduce the main features of the observed phenomena in the two scenarios.

2.2. Mathematical Framework

For the development of the mathematical modelling explained in the next section, we remark that we neglect the third dimension, since we do not have laboratory measurements of the movement of cells in the vertical direction. For this reason we consider the microfluidic chambers as 2D objects. Nowadays, the 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 allow to gain a deeper understanding of the interactions between cancer cells and the immune system. More generally, mathematical models can describe a broad variety of biological phenomena, including cell dynamics and cancer [28,29,30,31,32].

The movement of bacteria under the effect of a chemical substance has been widely studied in the last few decades, and numerous mathematical models have been proposed. As shown in [33], 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: particles models, hybrid (multiscale) models and macroscopic models.

In the present paper, the population density is assumed to be as a whole, thus macroscopic models of partial differential equations are considered. In particular, 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: a 1D KS-like model and a 1D GA-like model. The modelling here applied is described in the next Section 2.2.1.

2.2.1. The Model

Here, we introduce a mathematical model that aims at describing the behavior of two populations of cells coexisting together: tumoral cells T and immune cells (macrophages) M. We underline that the setting here considered can be made 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, interactions with chemoattractants, interactions 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 [17], ad hoc transmission conditions were introduced to ensure 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 rgw observed data.

Cancer cells T produce chemical signals, called $\phi$, activating the immune response of M and influencing their behavior. Moreover, we take into account the presence of cytokines $\omega$ (produced by M), acting as a chemical killer of cancer cells. Mainly referring to the KS model, the model here considered in the 2D chambers reads as:

$$\left\{\begin{array}{c}\frac{\partial }{\partial t}T=D_T\Delta T-{\lambda }_T\left(\omega \right)T,\hfill \\ \frac{\partial }{\partial t}M=D_M\Delta M-div(\chi \left(\phi \right)M\nabla \phi ),\hfill \\ \\ \frac{\partial }{\partial t}\phi =D_{\phi }\Delta \phi +{\alpha }_{\varphi }T-{\beta }_{\phi }\phi ,\hfill \\ \\ \frac{\partial }{\partial t}\omega =D_{\omega }\Delta \omega +{\alpha }_{\omega }M-{\beta }_{\omega }\omega .\hfill \\ \end{array}\right.$$

(3)

The system above describes the dynamics of the two cell species and the diffusion of the chemoattractant, and it needs to be complemented with suitable initial conditions and boundary conditions for the cells and the chemoattractant concentrations, as will be specified in the following paragraphs.

In particular, bearing in mind the first scenario (treated tumor), the system above describes the following situation: tumor cells T produce a chemical substance $\phi$ attracting immune cells M, that start migrating towards them. In this case, the tumor does not proliferate, since it is treated by a chemotherapy medicine, thus we do not include this feature in the model. Immune cells do not seem to proliferate during the experiment, thus we neglect this feature in the model. Note that, although here we are neglecting the proliferation of tumor, a tumor growth is observed experimentally in the second scenario (untreated tumor). This feature can be easily added to the model by putting a linear or logistic source term in the equation for the tumor.

Immune cells also produce a chemical substance $\omega$ which should be responsible for tumor killing. Therefore, in the first equation of the system (3), besides the diffusion term we can find $-{\lambda }_T\left(\omega \right)T$, representing the tumor suppression operated by immune cells. We underline that at the moment we have no information from the biologists about the real killing rate in the microchip environment induced by $\omega$, but we decided to introduce it in order to include this phenomenon in the model qualitatively.

In the second equation, in addition to the diffusion term we have the chemotactic term $f=\chi \left(\phi \right)\nabla \phi$ due to the presence of the chemical substance $\phi$ produced by the tumor.

In order to define the action of the cytotoxin $\omega$ produced by immune cells (which determines the death of cancer cells), we introduce the function ${\lambda }_T\left(\omega \right)$:

$${\lambda }_T\left(\omega \right)=\frac{S\omega }{\gamma +\omega },$$

(4)

where S is the maximum secretion rate of the cytotoxin by the immune cells and $\gamma$ the equivalent Michaelis constant associated with the production, as described in [33]. A lot of effort has been devoted to finding a biologically accurate expression for the chemotaxis function $\chi \left(\phi \right)$ representing the chemotactic sensitivity of immune cells; here, we refer the form suggested in [34] by Lapidis and Schiller:

$$\chi \left(\phi \right)=\frac{k_1}{{(k_2+\phi )}^2},$$

(5)

where $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. Note that we mainly refer to [33] for the values of the parameters $k_1$, $k_2$, and all the parameters of the model 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, with the label c indicating the channels. Then, we consider two possible models in the channels: a diffusive one:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}\frac{\partial }{\partial t}{T}_c=D_T{\partial }_{xx}{T}_c-{\lambda }_T\left(\omega \right)T_c,\hfill \\ \\ \frac{\partial }{\partial t}{M}_c=D_M{\partial }_{xx}{M}_c-{\partial }_x\left(M_c{f}_c\right),\hfill \\ \\ \frac{\partial }{\partial t}{\phi }_c=D_{\phi }{\partial }_{xx}{\phi }_c+{\alpha }_{\phi }{T}_c-{\beta }_{\phi }{\phi }_c,\hfill \\ \\ \frac{\partial }{\partial t}{\omega }_c=D_{\omega }{\partial }_{xx}{\omega }_c+{\alpha }_{\omega }{M}_c-{\beta }_{\omega }{\omega }_c,\hfill \\ \end{array}\right.\hfill \end{array}$$

(6)

and a hyperbolic one:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\partial }_t{T}_c+{\partial }_x{v}_c^T=-{\lambda }_T\left(\omega \right)T_c,\hfill \\ \\ {\partial }_t{v}_c^T+D_T{\partial }_x{T}_c=-v_c^T,\hfill \\ \\ {\partial }_t{\omega }_c=D_{\omega }{\partial }_{xx}{\omega }_c+{\alpha }_{\omega }{T}_c-{\beta }_c{\omega }_c,\hfill \\ \\ {\partial }_x{M}_c+{\partial }_t{v}_c^M=0,\hfill \\ \\ {\partial }_t{v}_c^M+D_M{\partial }_x{M}_c=M_c{f}_c-v_c^M,\hfill \\ \\ {\partial }_t{\phi }_c=D_{\phi }{\partial }_{xx}{\phi }_c+{\alpha }_{\phi }{T}_c-{\beta }_{\phi }{\phi }_c,\hfill \end{array}\right.\hfill \end{array}$$

(7)

with $T_c$ and $M_c$, respectively, as the density of tumor and immune cells; $f_c=\chi \left({\phi }_c\right){\partial }_x{\phi }_c$ and $v_c^T$ and $v_c^M$, respectively, as the average flux of tumor cells $T_c$ and immune cells $M_c$ in the channels. Note that the 1D-doubly parabolic model (6) is the one-dimensional version of the system (3), while (7) is the 1D-hyperbolic-parabolic model inspired by GA model [5]. We remark that, for the hyperbolic-parabolic system (7), we also need to assign initial and boundary conditions for the flux v.

2.2.2. The Simplified Model

To present the numerical scheme, we write a simpler model with respect to (3) but share the main features of it:

$$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\partial }_{t}u\hfill & =& D_u\Delta u-\mathrm{div}\mathbf{F}+g(x,y,t,u),\hfill \\ {\partial }_t\varphi \hfill & =& D_{\varphi }\Delta \varphi +au-b\varphi ,\hfill \end{array}\right.\end{array}$$

(8)

with u as the density of individuals, $\varphi$ as the density of chemoattractant, and with $\mathbf{F}=\mathbf{u}\mathbf{f}$. From now on, the two components of the drift term $\mathbf{f}=\chi \left(\varphi \right)\nabla \varphi$ will be indicated as:

$$\mathbf{f}\left(\mathbf{x},\mathbf{y},\mathbf{t}\right):=\left(\begin{array}{c}{f}^x\left(x,y,t\right)\\ f^y\left(x,y,t\right)\end{array}\right).$$

For the mono-dimensional channel, in order to make the explanation of the numerical approximation easier, we consider simpler models sharing the same characteristics of the models in (6) and (7), which read as:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\partial }_t{u}_c=D_{u_c}{\partial }_{xx}u-{\partial }_x{F}_c+g(x,t,u),\hfill \\ {\partial }_t{\varphi }_c=D_{{\varphi }_c}{\partial }_{xx}{\varphi }_c+a_c{u}_c-b_c{\varphi }_c\hfill \end{array}\right.\hfill \end{array}$$

(9)

or

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\partial }_t{u}_c+{\partial }_x{v}_c=g(x,t,u_c),\hfill \\ {\partial }_t{v}_c+{\lambda }_c^2{\partial }_x{u}_c=F_c-v_c,\hfill \\ {\partial }_t{\varphi }_c=D_{{\varphi }_c}{\partial }_{xx}{\varphi }_c+a_c{u}_c-b_c{\varphi }_c,\hfill \end{array}\right.\hfill \end{array}$$

(10)

with $F_c=u_c{f}_c$. The systems above have to be complemented with smooth initial conditions for the unknowns $u,\varphi$ and also v for system (10); initial data will be specified in Section 3.3. On the boundary, we consider for all the quantities homogeneous Neumann conditions, so we assume no-flux boundary conditions.

Monotonicity conditions. We also mention at this point that model (10) requires an analytical monotonicity criteria; see [35]. For a linear convection term $A{u}_c$ and a linear source term $g=B{u}_c$, the sub-characteristic criteria

$$\left|\frac{A}{{\lambda }_c}\right|-B\le 1,$$

must be satisfied in order for the quantity $u_c$ to be non-negative. Otherwise, having a negative $u_c$ would lead to unphysical solutions.

With regard to our former model (7), that would mean we have for the immune cell density M the monotonicity condition:

$$\frac{k_1}{{\left(k_2+{\phi }_c\right)}^2}\left|{\partial }_x{\phi }_c\right|\le \sqrt{D_M}.$$

(11)

This needs 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. Then, we aim at extending the no-flux boundary conditions to more general ones.

In the following section, in order to discuss mass-conservation we restrict our study on a 2D closed rectangular domain named ${\Omega }_{nc}$ (i.e., box without the channels) and on single lines not connected at the outer boundaries, as indicated by $C_{nc}$.

2.3. Outer Boundary and Interface Conditions for the Models with Null Source Term G

From now on, we consider a further simplified version of the models presented above putting the source term g equal to zero to discuss mass conservation.

2.3.1. Boundary Conditions for the 2d Doubly Parabolic Model (8) with $G=0$.

Here, we consider the mass conservation of cell density u for zero-flux conditions at the outer boundaries of a rectangular closed domain ${\Omega }_{nc}$ for the 2D model (8) with a null source term.

Indeed, since our model describes the migration of cells by both diffusion and chemoattractant effects, physically speaking the mass of cells must be preserved in the absence of the creation and destruction of cells.

For this reason, we assume a no-flux condition for the density u and the chemoattractant $\varphi$. Since we define the source term $\mathbf{f}$ as a product function of $▿\varphi$, we get the no-flux conditions:

$$\begin{array}{cc}\hfill \mathbf{F}\left(\mathbf{x},\mathbf{y},\mathbf{t}\right){\mathit{n}|}_{\partial {\Omega }_{nc}}& =0,\mathbf{▿}\mathit{u}\left(\mathit{x},\mathit{y},\mathit{t}\right){\mathit{n}|}_{\delta {\Omega }_{nc}}=0,(x,y)\in \partial {\Omega }_{nc},\end{array}$$

(12)

These guarantee mass conservation.

2.3.2. Interface between 2D-1D Models in (8) and (9)

We remark that, for simplicity, the following description focuses on the left part of the domain—i.e., ${\Omega }_l$ and its connection with a single microchannel represented by the interval $[0,L]$. However, the numerical treatment holds analogously on the right side of the domain and also to multiple channels.

In the 2D left box ${\Omega }_l$, the position of node $1L$ is at $x=L_y,y\in [a,b]$ and for the 1D domain, represented by a given number of channels C, node $1L$ is placed at $x=0$ (left endpoint of the channel); see Figure 2.

In order to ensure the conservation of the total mass when the mass-exchange occurs, we introduce suitable transmission conditions at the interface between 2D and 1D domains. Then, we consider the simplified models (9) or (10) for 1D channels, with $g=0$.

In particular, we have to prescribe the flux conservation at the 2D-1D interfaces, since we cannot lose or gain any cells during the passage through a node.

The conservation condition for u reads as:

$${\displaystyle \frac{d}{dt}{\int }_{{\Omega }_l}u\left(x,y,t\right)d{\Omega }_l+\frac{d}{dt}{\int }_0^L{u}_c\left(x,t\right)dx=0,}$$

and it rewrites as:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle 0={\oint }_{\partial {\Omega }_l}\left(D_u\mathbf{▿}\mathit{u}\left(\mathit{x},\mathit{y},\mathit{t}\right)-\mathit{F}\left(\mathit{x},\mathit{y},\mathit{t}\right)\right)\mathit{n}dS}\hfill \\ & +(D_{u_c}{\partial }_x{u}_c-F_c(x,t)){|}_L^0\hfill \end{array}\end{array}$$

by using the divergence theorem in the first integral. With our analytical boundary conditions, the integral vanishes, except at the boundary where the node is positioned.

We remark that attention has to be paid to $\mathit{n}$ being the outer normal of the domain.

Thanks to the boundary conditions (12), we get the condition:

$${\displaystyle {\int }_a^b\left(D_u{\partial }_{x}u\left(L_x,y,t\right)-u(L_x,y,t)f^x(L_x,y,t\right))dy=D_{u_c}{\partial }_x{u}_c\left(0,t\right)-F_c\left(0,t\right).}$$

(13)

Note that in this case, the evaluation of the integrand at the right endpoint L is discarded, since we are only considering the junction connecting ${\Omega }_l$ with the left endpoint $x=0$ of the microchannel (see Figure 2).

Now, we impose Kedem–Katchalsky (KK) [20] conditions describing the conservation of the flux through a node (see also [21] for the numerical treatment of these conditions). In particular, at the interface between the left chamber and the channels, we have (on the left of node $1L$ in Figure 2):

$$D_u{\partial }_{x}u\left(L_x,y,t\right)-u\left(L_x,y,t\right)f^x\left(L_x,y,t\right)=K\left(u_c\left(0,t\right)-u(L_x,y,t)\right)\mathrm{for}y\in \left[a,b\right]$$

(14)

On the right of node $1L$, we have:

$${\displaystyle D_{u_c}{\partial }_x{u}_c\left(0,t\right)-F_c\left(0,t\right)=K{u}_c(0,t)(b-a)-{\int }_a^{b}u(L_x,y,t)dy.}$$

(15)

Thanks to conditions (14) and (15), we are guaranteed to have the flux conservation (13), 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 3 and in Section 3.1.2.

2.3.3. Interface between 2D-1D Models in (8)–(10)

In this section, we describe the combination of 2D parabolic-1D hyperbolic model in order to describe the dynamics with a hyperbolic model in the channels. Further care has to be given in order to keep some important properties which ensure the consistency and non-negativity of numerical solutions when connecting both models.

Now, the transmission conditions for the switch from ${\Omega }_l$ to $C=[0,L]$ are derived in this case. For the mass conservation of u, we impose the condition:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle 0=\frac{d}{dt}{\int }_{{\Omega }_l}u\left(x,y,t\right)d{\Omega }_l+\frac{d}{dt}{\int }_0^L{u}_c\left(x,t\right)dx}\hfill \\ & {\displaystyle ={\int }_{{\Omega }_l}\left(D_u▵u\left(x,y,t\right)-\mathrm{div}\mathit{F}\left(x,y,t\right)\right)d{\Omega }_l+{\int }_0^L-{\partial }_{x}v\left(x,t\right)dx}\hfill \\ ⟹\hfill & {\displaystyle {\oint }_{\partial {\Omega }_l}\left(D_u\mathbf{▿}\mathit{u}\left(\mathit{x},\mathit{y},\mathit{t}\right)-\mathit{F}(\mathit{x},\mathit{y},\mathit{t})\mathit{n}\right)dS+v\left(0,t\right)=0.}\hfill \end{array}\end{array}$$

Note that in the above formula, we have $v\left(L,t\right)=0$ because we are looking at left interface (node $1L$). Then, we finally get:

$${\displaystyle {\int }_{a_1}^{b_1}\left(D_u{\partial }_{x}u\left(L_x,y,t\right)-u\left(L_x,y,t\left)f^x\right(L_x,y,t\right)\right)dy=-v\left(0,t\right).}$$

(16)

Now, we impose the KK-condition at the interface:

$$D_u{\partial }_{x}u\left(L_x,y,t\right)-u(L_x,y,t)f^x\left(L_x,y,t\right)=K\left(u_c\left(0,t\right)-u\left(L_x,y,t\right)\right)\mathrm{for}y\in \left[a,b\right]$$

Then, (16) reads as:

$$\begin{array}{ccc}\hfill v\left(0,t\right)& =& {\displaystyle -K(b-a)u_c\left(0,t\right)+K{\int }_a^{b}u\left(L_x,y,t\right)dy.}\hfill \end{array}$$

(17)

3. Numerical Approximation

In this section, we describe the numerical approximation of the adopted models: 2D-doubly parabolic, 1D-doubly parabolicm and 1D-hyperbolic-parabolic. We define equispaced $x_i:=i▵x$, $y_j:=j▵y$, and $t_n:=n▵t$ with $▵x$, $▵y$, $▵t>0$, and $i=0,\dots ,N_x+1$, $j=0,\dots ,N_y+1$; for channel $[0,L]$, we discretize it as $x_i=i▵x$, with $i=0,\dots ,N$. For a more structured presentation, we introduce the operators:

$$\begin{array}{c}{\delta }_x^2{u}_{i,j}^n:=u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n,{\delta }_y^2{u}_{i,j}^n:=u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n,\hfill \\ {\delta }_x^0{u}_{i,j}^n:=u_{i+1,j}^n-u_{i-1,j}^n,{\delta }_y^0{u}_{i,j}^n:=u_{i,j+1}^n-u_{i,j-1}^n,\hfill \\ {\delta }_x^1{u}_{i,j}^n:=u_{i+1,j}^n-u_{i,j}^n,{\delta }_y^1{u}_{i,j}^n:=u_{i,j+1}^n-u_{i,j}^n.\hfill \end{array}$$

Remark 2.

Note that special attention has to be paid also to the stiffness induced by the source term $g\left(x,y,t,u\right)$. To overcome this issue, implicit methods can be used, such as the Crank–Nicolson method, which is associated with a $\Delta t$ that is small enough.

Mass-preserving and (numerically)-positivity-preserving approximations will be developed in the present section. In the following, we will neglect the label c to make the reading easier and we will make distinctions only when necessary.

3.1. The Parabolic-Parabolic Case

Here, we introduce a numerical scheme for the doubly parabolic systems (8) and (9) for $g=0$.

For the discretization of equations in a 2D system (8) in the interior points of the domain—i.e., for $i=1,\dots ,N_x,j=1,\dots ,N_y$, we define a finite difference discretization both for u and $\varphi$:

$$\begin{array}{ccc}{u}_{i,j}^{n+1}\hfill & =& u_{i,j}^n+D_u\frac{▵t}{2}\left[\frac{{\delta }_x^2(u_{i,j}^n+u_{i,j}^{n+1})}{▵x^2}+\frac{{\delta }_y^2(u_{i,j}^n+u_{i,j}^{n+1})}{▵y^2}\right]\hfill \\ \hfill & & -▵t({\delta }_x^0{F}_{i,j}^{x,n}+{\delta }_y^0{F}_{i,j}^{y,n}),\hfill \end{array}$$

(18)

$$\begin{array}{ccc}{\varphi }_{i,j}^{n+1}\hfill & =& {\varphi }_{i,j}^n+D_{\varphi }\frac{▵t}{2}\left[\frac{{\delta }_x^2\left({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1}\right)}{▵x^2}+\frac{{\delta }_y^2\left({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1}\right)}{▵y^2}\right]\hfill \\ \hfill & & +\frac{▵t}{2}a(u_{i,j}^n+u_{i,j}^{n+1})-\frac{▵t}{2}b({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1}).\hfill \end{array}$$

(19)

Note that $F_i^{x,n}=\chi ({\varphi }_{i,j}^n)u_{i,j}^n{\delta }_x^0{\varphi }_{i,j}^n$ with ${\delta }_x^0{\varphi }_{i,j}^n$ a centered second order approximation of ${\varphi }_x$.

For a 1D system (9), in the interior points of the channel C we apply the Crank–Nicolson scheme, as above:

$$\begin{array}{ccc}{u}_i^{c,n+1}\hfill & =& u_i^{c,n}+D_{u_c}\frac{\Delta t}{2}\frac{{\delta }_x^2\left(u_i^{c,n}+u_i^{c,n+1}\right)}{▵x^2}-\Delta t{\delta }_x^0{F}_{c,i}^n,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}{\varphi }_i^{c,n+1}\hfill & =& {\varphi }_i^{c,n}+D_{{\varphi }_c}\frac{▵t}{2}\frac{{\delta }_x^2\left({\varphi }_i^{c,n}+{\varphi }_i^{c,n+1}\right)}{▵x^2}\hfill \\ \hfill & & +\frac{▵t}{2}{a}_c(u_i^{c,n}+u_i^{c,n+1})-\frac{▵t}{2}{b}_c({\varphi }_i^{c,n}+{\varphi }_i^{c,n+1}).\hfill \end{array}$$

(21)

Remark 3.

We remark that here, for simplicity, upwinding terms needed to avoid meshgrid restrictions caused by a large Péclet number—see, for instance [36]—are not included in the schemes (18) and (20). In this case, assuming a small $\Delta x$ and $\Delta y$ will be enough to avoid oscillations.

The upwinding be added in the implemented scheme described in Section 3.3.

Stability condition. By using the Von-Neumann stability analysis on the linearized problem, we derive the following condition for the Crank–Nicolson schemes above, [37]:

$$D_u\frac{▵t}{▵x^2}\le 1\mathrm{for}1\mathrm{D},D_u\frac{▵t}{▵x^2}+D_u\frac{▵t}{▵y^2}\le 1\mathrm{for}2\mathrm{D}.$$

(22)

In the following, we present the discretization of the boundary and transmission conditions to complete the numerical schemes.

3.1.1. Discretization of the Outer Boundary Conditions for the Doubly Parabolic Problem

Since a qualitative characteristic of the model is the preservation of total mass for zero-flux boundary conditions, the first step is to ensure the mass preservation at each time step in a closed 1D line $C_{nc}$ and, analogously, in a closed 2D chamber, namely ${\Omega }_{nc}$. 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 following to complete the equations.

Boundary conditions for the density of individuals u in 1D model (9).

Here, we consider Neumann boundary conditions $\frac{\partial u}{\partial x}(x,t)=0$. If we discretized it with a (standard) forward finite difference scheme:

$${\partial }_{x}u\left(0\right)=\frac{u_1^{c,n}-u_0^{c,n}}{\Delta x},$$

(23)

the mass will not be preserved over time, as shown in the numerical Example 1. Therefore, in order to have mass preservation, we use the central finite difference approximation with a ghost cell:

$${\partial }_{x}u\left(0\right)=\frac{u_1^{c,n}-u_{-1}^{c,n}}{2\Delta x}.$$

(24)

At the outer boundaries at the first ($x_0=0$) and last ($x_{N+1}=L$) endpoint of the channels, we assign the numerical schemes:

$$u_0^{c,n+1}=u_0^{c,n}+D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_0^{c,n}+u_0^{c,n+1}\right)-\frac{▵t}{▵x}(F_{c,1}^n+F_{c,0}^n)$$

(25)

and

$$u_{N+1}^{c,n+1}=u_{N+1}^{c,n}-D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_N^{c,n}+u_N^{c,n+1}\right)+\frac{▵t}{▵x}(F_{c,N}^n+F_{c,N+1}^n),$$

(26)

with $F_{c,0}^n=0$ in (25) and $F_{c,N+1}^n=0$ in (26), since we are imposing zero flux at the boundaries. The boundary conditions above are mass-preserving by construction, as stated in the following proposition.

Proposition 1.

The 1D numerical scheme at the internal points introduced above, namely:

$$u_i^{c,n+1}=u_i^{c,n}+D_{u_c}\frac{▵t}{2▵x^2}{\delta }_x^2(u_i^{c,n}+u_i^{c,n+1})-\frac{▵t}{2▵x}{\delta }_x^0{F}_{c,i}^n,$$

(27)

endowed with boundary conditions (25) and (26), is mass-preserving by construction, since it is obtained by imposing ${\mathcal{I}}^{n+1}-{\mathcal{I}}^n=0$.

Proof. 

The mass conservation over time on the closed domain ${\Omega }_{nc}$ reads as:

$$\begin{array}{c}\hfill I\left(t\right)={\int }_{{\Omega }_{nc}}u\left(t,x,y\right)d{\Omega }_{nc}={\int }_{{\Omega }_{nc}}u\left(0,x,y\right)d{\Omega }_{nc}=I\left(0\right).\end{array}$$

(28)

Now, applying a quadrature rule for the numerical integration:

$$\begin{array}{c}\hfill {\mathcal{I}}^n\approx {\int }_{{\Omega }_{nc}}u\left(t,x,y\right)d{\Omega }_{nc}.\end{array}$$

(29)

We need to ensure that ${\mathcal{I}}^{n+1}={\mathcal{I}}^n$.

For the numerical integration, different quadrature formulas can be used; in particular, we applied closed Newton–Cotes formulas to take into account the values at the boundaries. By using the trapezoidal rule with an integration error of $\mathcal{O}\left(▵x^2\right)$ and imposing the equality ${\mathcal{I}}^{n+1}-{\mathcal{I}}^n=0$ in the 1D case, one obtains:

$$\begin{array}{c}\hfill ▵x\left({\displaystyle \frac{u_0^{c,n+1}}{2}-\frac{u_0^{c,n}}{2}+\sum _{i=1}^N\left(u_i^{c,n+1}-u_i^{c,n}\right)+\frac{u_{N+1}^{c,n+1}}{2}-\frac{u_{N+1}^{c,n}}{2}}\right)=0.\end{array}$$

Using the numerical scheme (20) for $u_i^{c,n+1}$ in (30) for $i=1,\dots ,N$ we have:

$$\begin{array}{c}\hfill \begin{array}{ccc}▵x\hfill & \left(\frac{u_0^{c,n+1}-u_0^{c,n}}{2}\right.\hfill & +\frac{D_{u_c}▵t}{2▵x^2}\sum _{i=1}^N{\delta }_x^2\left(u_i^{c,n}+u_i^{c,n+1}\right)\hfill \\ & & \left.{\displaystyle -\frac{▵t}{2▵x}\sum _{i=1}^N{\delta }_x^0{F}_{c,i}^n+\frac{u_{N+1}^{c,n+1}-u_{N+1}^{c,n}}{2}}\right)=0.\hfill \end{array}\end{array}$$

Now, using the definition of ${\delta }_x^0{F}_{c,i}^n$ in the sum, the above formula becomes:

$$\begin{array}{c}\hfill \begin{array}{c}{u}_0^{c,n+1}-u_0^{c,n}-D_{u_c}\frac{▵t}{▵x^2}\left(u_1^{c,n}+u_1^{c,n+1}-(u_0^{c,n}+u_0^{c,n+1})\right)+\frac{▵t}{▵x}\left(F_{c,0}^n+F_{c,1}^n\right)\hfill \\ +u_{N+1}^{c,n+1}-u_{N+1}^{c,n}-D_{u_c}\frac{▵t}{▵x^2}\left(u_N^{c,n}+u_N^{c,n+1}-(u_{N+1}^{c,n}+u_{N+1}^{c,n+1})\right)\hfill \\ -\frac{▵t}{▵x}\left(F_{c,N}^n+F_{c,N+1}^n\right)=0.\hfill \end{array}\end{array}$$

We can now compute the values for both $u_0^{n+1}$ and $u_{N+1}^{n+1}$ so that the term equals zero. By collecting values from nearby stencils together (otherwise we obtain an error of $\mathcal{O}(▵x)$, which can be verified by Taylor expansion), at the outer boundaries of 1D domain we obtain the schemes (25) and (26). Note that in the formula above, by imposing homogeneous Neumann boundary conditions we have $f_0^n=0$ and $f_{N+1}^n=0$.

We also remark that besides imposing the stability condition (22), $\Delta x$ has to be small enough in order to ensure the positivity of the scheme for taking into account the possibility of having a negative term brought by $f_1^n$ or $f_N^n$ in (25) and (26), respectively. □

Remark 4.

Although here we are not proving that, the scheme (27) is in practice second order in the space up to the boundaries, since it can be equivalently obtained using the second-order approximation of the first derivative including a ghost cell reported in (24).

We underline that for $\mathbf{f}\ne 0$, even using formula (24), the approach with the discrete integral equation is necessary for ensuring the mass preservation. Furthermore, by using a different numerical integration scheme, we can achieve different mass-preserving boundary conditions of a higher order.

Boundary conditions for the density of individuals u in 2D model (8).

Let us assume there are no-flux conditions at the boundaries. Then, in this case we consider the 2D closed domain ${\Omega }_{nc}$.

Using the mass-preserving property argument, we compute boundary conditions for the corners and the top and bottom boundaries of ${\Omega }_{nc}$. By applying them with the numerical method (18) into ${\mathcal{I}}^{n+1}-{\mathcal{I}}^n=0$, we obtain the expression:

$$\frac{▵t▵x}{4}\left({\displaystyle -4▵t\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\left({\delta }_x^0{F}_{i,j}^x+{\delta }_y^0{F}_{i,j}^y\right)}\right)=0,$$

since the terms in u cancel. By plugging in it the expression of the central in rgw space second-order finite difference scheme ${\delta }_x^0{f}_j^{x,n}$ for $\mathrm{div}\left(\mathbf{f}\right)$, we obtain:

$$\begin{array}{cc}& {\displaystyle \frac{1}{▵y}\sum _{i=1}^{N_x}\left(F_{i,N_y+1}^{y,n}+F_{i,N_y}^{y,n}-F_{i,1}^{y,n}-F_{i,0}^{y,n}\right)}\\ & {\displaystyle +\frac{1}{▵x}\sum _{j=1}^{N_y}\left(F_{N_x+1,j}^{x,n}+F_{N_x,j}^{x,n}-F_{1,j}^{x,n}-F_{0,j}^{x,n}\right)=0.}\end{array}$$

Now, we can distribute the remaining values to the boundaries in the same way as above for the 1D-parabolic case. Therefore, we obtain the following mass-preserving boundary conditions for the corners:

$$\left\{\begin{array}{ccc}{u}_{0,0}^{n+1}\hfill & =\hfill & u_{0,0}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^1\left(u_{0,0}^n+u_{0,0}^{n+1}\right)+D_u\frac{▵t}{2▵y^2}{\delta }_y^1\left(u_{0,0}^n+u_{0,0}^{n+1}\right),\hfill \\ \\ u_{N_x+1,0}^{n+1}\hfill & =\hfill & u_{N_x+1,0}^n-D_u\frac{▵t}{2▵x^2}{\delta }_x^1\left(u_{N_x,0}^n+u_{N_x,0}^{n+1}\right)\hfill \\ & \hfill & +D_u\frac{▵t}{2▵y^2}{\delta }_y^1\left(u_{N_x+1,0}^n+u_{N_x+1,0}^{n+1}\right),\hfill \\ u_{0,N_y+1}^{n+1}\hfill & =\hfill & u_{0,N_y+1}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^1\left(u_{0,N_y+1}^n+u_{0,N_y+1}^{n+1}\right)\hfill \\ & \hfill & -D_u\frac{▵t}{2▵y^2}{\delta }_y^1\left(u_{0,N_y}^n+u_{0,N_y}^{n+1}\right),\hfill \\ u_{N_x+1,N_y+1}^{n+1}\hfill & =\hfill & u_{N_x+1,N_y+1}^n-D_u\frac{▵t}{2▵x^2}{\delta }_x^1\left(u_{N_x,N_y+1}^n+u_{N_x,N_y+1}^{n+1}\right)\hfill \\ & \hfill & -D_u\frac{▵t}{2▵y^2}{\delta }_y^1\left(u_{N_x+1,N_y}^n+u_{N_x+1,N_y}^{n+1}\right).\hfill \end{array}\right.$$

(30)

For the edges of the box ${\Omega }_{nc}$, $i=1,\dots ,N_x$ and $j=1,\dots ,N_y$, we have:

$$\left\{\begin{array}{ccc}{u}_{i,0}^{n+1}\hfill & =\hfill & u_{i,0}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^2\left(u_{i,0}^n+u_{i,0}^{n+1}\right)+D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{i,0}^n+u_{i,0}^{n+1}\right)\hfill \\ \hfill & \hfill & -\frac{▵t}{▵y}{F}_{i,1}^{y,n},\hfill \\ u_{i,N_y+1}^{n+1}\hfill & =\hfill & u_{i,N_y+1}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^2\left(u_{i,N_y+1}^n+u_{i,N_y+1}^{n+1}\right)\hfill \\ & \hfill & -D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{i,N_y}^n+u_{i,N_y}^{n+1}\right)+\frac{▵t}{▵y}{F}_{i,N_y}^{y,n},\hfill \\ \\ u_{0,j}^{n+1}\hfill & =\hfill & u_{0,j}^n+D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{0,j}^n+u_{0,j}^{n+1}\right)+D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{0,j}^n+u_{0,j}^{n+1}\right)\hfill \\ \hfill & \hfill & -\frac{▵t}{▵x}{F}_{1,j}^{x,n},\hfill \\ u_{N_x+1,j}^{n+1}\hfill & =\hfill & u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,j}^n+u_{N_x,j}^{n+1}\right)\hfill \\ & \hfill & +D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)+\frac{▵t}{▵x}{F}_{i,N_y}^{y,n}.\hfill \end{array}\right.$$

(31)

We underline that in the above formulas the terms $F^{x,n}$ and $F^{y,n}$ on the edges of the domain cancel because of the homogeneous Neumann boundary condition on ${\varphi }_x$ and ${\varphi }_y$.

Boundary conditions for the density of chemoattractant $\varphi$ in 1D model (9).

For the computation of the conditions at the outer boundaries for the chemoattractant ${\varphi }_c$ in the 1D-doubly parabolic model we proceed as above, but neglect the source term $a_c{u}_c-b_c{\varphi }_c$ to obtain boundary conditions that are mass-preserving.

Proceeding as above, we achieve the following mass-preserving boundary conditions for the chemoattractant:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\varphi }_0^{c,n+1}\hfill & ={\varphi }_0^{c,n}+D_{{\varphi }_c}\frac{▵t}{▵x^2}\left({\delta }_x^1{\varphi }_0^{c,n}+{\delta }_x^1{\varphi }_0^{c,n+1}\right)\\ {\varphi }_{N+1}^{c,n+1}\hfill & ={\varphi }_{N+1}^{c,n}-D_{{\varphi }_c}\frac{▵t}{▵x^2}\left({\delta }_x^1{\varphi }_N^{c,n}+{\delta }_x^1{\varphi }_N^{c,n+1}\right)\end{array}\right.\end{array}$$

(32)

This, under the CFL condition (22), is second-order accurate and positivity-preserving up to the boundaries. The parabolic equation in the interior points is solved using an implicit-explicit method:

$${\varphi }_i^{c,n+1}={\varphi }_i^{c,n}+D_{{\varphi }_c}\frac{▵t}{2▵x^2}{\delta }_x^2\left({\varphi }_i^{c,n}+{\varphi }_i^{c,n+1}\right).$$

(33)

Boundary conditions for the density of chemoattractant $\varphi$ in 2D model (8).

Reasoning as above, by applying the numerical method (19) we obtain the following boundary condition for the chemoattractant at the corners:

$$\left\{\begin{array}{ccc}{\varphi }_{0,0}^{n+1}\hfill & =\hfill & {\varphi }_{0,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{0,0}^n+{\varphi }_{0,0}^{n+1})+2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{0,0}^n+{\varphi }_{0,0}^{n+1})\hfill \\ \\ {\varphi }_{N_x+1,0}^{n+1}\hfill & =\hfill & {\varphi }_{N_x+1,0}^n-D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x+1,0}^n+{\varphi }_{N_x+1,0}^{n+1})\hfill \\ \hfill & \hfill & +D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{N_x+1,0}^n+{\varphi }_{N_x+1,0}^n)\hfill \\ \\ {\varphi }_{N_x+1,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{Nx+1,N_y+1}^n-D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x,N_y+1}^n+{\varphi }_{N_x,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -D_{\varphi }\frac{▵t}{▵y^2}{\delta }_x^1({\varphi }_{N_x+1,N_y+1}^n+{\varphi }_{N_x+1,N_y}^n)\hfill \\ \\ {\varphi }_{0,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{0,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{1,N_y+1}^n+{\varphi }_{1,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{0,N_y}^n+{\varphi }_{0,N_y}^{n+1})\hfill \end{array}\right.$$

(34)

and

$$\left\{\begin{array}{ccc}{\varphi }_{i,0}^{n+1}\hfill & =\hfill & {\varphi }_{i,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^2({\varphi }_{i,0}^n+{\varphi }_{i,0}^{n+1})+2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{i,0}^n+{\varphi }_{i,0}^{n+1})\hfill \\ \\ {\varphi }_{i,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{i,N_y+1}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^2({\varphi }_{i,N_y+1}^n+{\varphi }_{i,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{i,N_y}^n+{\varphi }_{i,N_y}^{n+1})\hfill \\ {\varphi }_{0,j}^{n+1}\hfill & =\hfill & {\varphi }_{0,j}^n+D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^2({\varphi }_{0,j}^n+{\varphi }_{0,j}^{n+1})+2D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{0,j}^n+{\varphi }_{0,j}^{n+1})\hfill \\ \\ {\varphi }_{N_x+1,j}^{n+1}\hfill & =\hfill & {\varphi }_{N_x+1,j}^n+D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^2({\varphi }_{N_x+1,j}^n+{\varphi }_{N_x+1,j}^{n+1}).\hfill \\ \hfill & \hfill & -2D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x,j}^n+{\varphi }_{N_x,j}^{n+1})\hfill \end{array}\right.$$

(35)

We now have a complete numerical method to solve the system (8) on ${\Omega }_{nc}$ and the 1D version of it (9) on $C_{nc}$.

Let us now turn to the complete domain depicted in Figure 2 in order to develop mass-preserving transmission conditions at the nodes of the network. We remark that for the sake of clarity, for the development of the numerical transmission conditions here we consider just the junction—indicated as node $1L$—connecting the left box ${\Omega }_l$ and a single channel parametrized as an interval $C=[a,b]$.

3.1.2. Discretization of the Transmission Conditions for the 2D-1D Doubly Parabolic Case

The choice of suitable transmission conditions is crucial, since it should reflect the qualitative attributes of the analytical model.

Here, we use ghost values (24) in order to obtain the numerical boundary conditions, since in such a way mass-preserving and positivity-preserving properties are ensured.

By using the central approximation formula for $div\mathit{f}$ in the condition (14) on the left of node $1L$, we have:

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill & D_u{\partial }_{x}u\left(L_x,y,t\right)-F^x\left(L_x,y,t\right)& =& K\left(u_c\left(0,t\right)-u\left(L_x,y,t\right)\right)\mathrm{for}y\in \left[a,b\right].\hfill \end{array}\end{array}$$

Then we have:

$$\begin{array}{c}\hfill \begin{array}{cccc}& D_u\frac{u_{N_x+2,j}^n-u_{N_x,j}^n}{2▵x}& =& K\left(u_0^{c,n}-u_{N_x+1,j}^n\right)+F_{N_x+1,j}^{x,n}\hfill \end{array}\end{array}$$

and we get:

$$\begin{array}{c}\hfill \begin{array}{ccc}{u}_{N_x+2,j}^n\hfill & =& u_{N_x,j}^n+K\frac{2▵x}{D}\left(u_0^{c,n}-u_{N_x+1,j}^n\right)+\frac{2▵x}{D_u}{F}_{N_x+1,j}^{x,n}\end{array}\end{array}$$

(36)

for $j=j_a,\dots ,j_b$.

Moreover, using the central approximation for ${\partial }_x{u}_c$ in (15),

we finally get the formula:

$${\displaystyle u_{-1}^{c,n}=u_1^{c,n}-K\frac{2▵x}{D_{u_c}}(b-a)u_0^{c,n}+\frac{2▵x}{D_{u_c}}\Delta y\sum _{j=j_a}^{j_b}K(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1})-\frac{2▵x}{D_{u_c}}{F}_{c,0}^n.}$$

(37)

We now use the Ansatz to apply the ghost values into the 1D (27) and 2D (18) numerical schemes without specific chemotactic approximation, and use the discrete integral equation to determine the chemotactic term discretization.

Since we need to conserve the mass in each domain, but also in both connected ones, the expanded discrete integral equation is needed to compute the total mass over both domains.

Plugging the ghost values (36) and (37), respectively, into the numerical schemes (18) and (27), we get the conditions at the interface (node $1L$):

$$\begin{array}{ccc}{u}_{N_x+1,j}^{n+1}\hfill & =& u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,j}^n+u_{N_x,j}^{n+1}\right)+2K\frac{▵t}{▵x}\left(u_0^{c,n}-u_{N_x+1,j}^n\right)\hfill \\ \hfill & & +D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)\hfill \\ \hfill & & +2\frac{▵t}{▵x}{F}_{N_x+1,j}^{x,n}-▵t{\delta }_x^0{F}_{N_x+1,j}^{x,n}-▵t{\delta }_y^0{F}_{N_x+1}^{y,n},\hfill \end{array}$$

(38)

and

$$\begin{array}{ccc}{u}_0^{c,n+1}\hfill & =& u_0^{c,n}+D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_0^{c,n}+u_0^{c,n+1}\right)-2K\frac{▵t}{▵x}(b-a)u_0^{c,n}\hfill \\ \hfill & & {\displaystyle +2K\frac{▵t▵y}{▵x}\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)}\hfill \\ \hfill & & -2\frac{▵t}{▵x}{F}_{c,0}^n-\frac{▵t}{2\Delta x}{\delta }_x^0{F}_{c,0}^n.\hfill \end{array}$$

(39)

In particular, the conservation of the discrete total mass reads as:

$${\mathcal{I}}_{1\mathrm{D}}^{n+1}+{\mathcal{I}}_{2\mathrm{D}}^{n+1}-{\mathcal{I}}_{1\mathrm{D}}^n-{\mathcal{I}}_{2\mathrm{D}}^n=0,$$

(40)

Now, applying the conditions (38) and (39) with the other boundary conditions (31) and (25), we get:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle ▵x(-K\frac{▵t}{▵x}(b-a)u_0^{c,n}+K\frac{▵t▵y}{▵x}\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)}\hfill \\ & -\frac{▵t}{▵x}{F}_{c,0}^n-\frac{▵t}{2}{\delta }_x^0{F}_{c,0}^n+\frac{▵t}{2▵x}\left(F_{c,0}^n+F_{c,1}^n\right))\hfill \\ +\hfill & \frac{▵x▵y}{4}({\displaystyle 2\sum _{j=j_a}^{j_b}(2K\frac{▵t}{▵x}\left(u_0^{c,n}-u_{N_x+1,j}^n\right)+2\frac{▵t}{▵x}{F}_{N_x+1,j}^{x,n}-▵t{\delta }_x^0{F}_{N_x+1,j}^{x,n}}\hfill \\ & -▵t{\delta }_y^0{F}_{N_x+1,j}^{y,n})\hfill \\ & {\displaystyle -\frac{2▵t}{▵x}\sum _{j=j_{a_1}}^{j_{b_1}}\left(F_{N_x+1,j}^{x,n}+F_{N_x,j}^{x,n}\right)-\frac{2▵t}{▵y}\sum _{j=j_a}^{j_b}\left(F_{N_x+1,j}^{y,n}+F_{N_x,j}^{y,n}\right))=0,}\hfill \end{array}\end{array}$$

We can then obtain the following transmission conditions:

$$\begin{array}{ccc}{u}_0^{c,n+1}\hfill & =& \underset{\mathrm{same}\mathrm{as}\mathrm{for}\mathrm{BC}\mathrm{without}\mathrm{transmission}\mathrm{condition}}{\underbrace{u_0^{c,n}+D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_0^{c,n}+u_0^{c,n+1}\right)-\frac{▵t}{▵x}(F_{c,1}^n+F_{c,0}^n)}}\hfill \\ \hfill & & {\displaystyle -2K\frac{▵t}{▵x}(b-a)u_0^{c,n}+2K\frac{▵t}{▵x}\sum _{j=j_a}^{j_b}{u}_{N_x+1,j}^n+u_{N_x+1,j}^{n+1},}\hfill \end{array}$$

(41)

$$\begin{array}{ccc}{u}_{N_x+1,j}^{n+1}\hfill & =& u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,j}^n+u_{N_x,j}^{n+1}\right)\hfill \\ & & +D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)\hfill \\ \hfill & & +\frac{▵t}{▵x}\left(F_{N_x+1,j}^{x,n}+F_{N_x,j}^{x,n}\right)-\frac{▵t}{2▵y}\left(F_{N_x+1,j+1}^{y,n}-F_{N_x+1,j-1}^{y,n}\right)\hfill \\ \hfill & & +\underset{\mathrm{additional}\mathrm{term}\mathrm{for}\mathrm{transmission}\mathrm{condition},}{\underbrace{2K\frac{▵t}{▵x}\left(u_0^{c,n}-u_{N_x+1,j}^n\right)}}\hfill \end{array}$$

(42)

for $j=j_a,\dots ,j_b$.

Proceeding analogously as above, this approach leads to mass-preserving and positivity-preserving transmission conditions for the chemoattractant $\varphi$ as well. In particular, we have at the first and last endpoint, respectively:

$$\begin{array}{c}\hfill \begin{array}{ccc}{\varphi }_0^{c,n+1}\hfill & =& {\varphi }_0^{c,n}+D_{{\varphi }_c}\frac{▵t}{▵x^2}\left({\delta }_x^1{\varphi }_0^{c,n}+{\delta }_x^1{\varphi }_0^{c,n+1}\right)+▵t{a}_c{u}_0^{c,n}-▵t{b}_c{\varphi }_0^{c,n}\hfill \\ & & {\displaystyle -2K\frac{▵t}{▵x}(b-a){\varphi }_0^{c,n}+2K\frac{▵t}{▵x}{\int }_a^b\varphi \left(L_x,y,t_n\right)dy}\hfill \end{array}\end{array}$$

and

$$\begin{array}{c}\hfill \begin{array}{ccc}{\varphi }_{N_x+1,j}^{n+1}\hfill & =& {\varphi }_{N_x+1,j}^n+D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^2({\varphi }_{N_x+1,j}^n+{\varphi }_{N_x+1,j}^{n+1})\hfill \\ & & -2D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x,j}^n+{\varphi }_{N_x,j}^{n+1})\hfill \\ & & +a▵t{u}_{N_x+1,j}^n-b▵t{\varphi }_{N_x+1,j}^n+2K\frac{▵t}{▵x}\left({\varphi }_0^{c,n}-{\varphi }_{N_x+1,j}^n\right).\hfill \end{array}\end{array}$$

(43)

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 by construction the mass conservation as the original PDE. Numerically, the scheme also ensures the positivity preserving property under the monotonicity conditions discussed in Section 3.3.1.

3.2. The Hyperbolic-Parabolic Case

The second-order AHO scheme on a line was introduced in [14] 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 of second order reads as:

$$\left\{\begin{array}{ccc}{u}_i^{c,n+1}\hfill & =& u_i^{c,n}+\lambda \frac{▵t}{2▵x}{\delta }_x^2{u}_i^{c,n}-\left(\frac{▵t}{2▵x}-\frac{▵t}{4\lambda }\right){\delta }_x^0{v}_i^{c,n}+\frac{▵t}{4\lambda }{\delta }_x^0{F}_{c,i}^n\hfill \\ & & +\frac{▵t}{4}\left(g_{i-1}^n+2g_i^n+g_{i+1}^n\right),\hfill \\ \\ v_i^{c,n+1}\hfill & =& v_i^{c,n}-{\lambda }^2\frac{▵t}{2▵x}{\delta }_x^0{u}_i^{c,n}+\left(\frac{\lambda ▵t}{2▵x}-\frac{▵t}{4}\right){\delta }_x^2{v}_i^{c,n}+\frac{▵t}{4}{\delta }_x^2{F}_{c,i}^n\hfill \\ \hfill & & +\lambda \frac{▵t}{4}\left(g_{i-1}^n-g_{i+1}^n\right),\hfill \end{array}\right.$$

(44)

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 [14]. For this reason we use the so-called AHO (Asymptotic Higher Order) schemes (see [18] 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 the third order of accuracy and converges towards a numerical solution with $v=0$, while preserving the mass.

3.2.1. Discretization of Transmission Conditions for the 2D-Doubly Parabolic and 1D-Hyperbolic-Parabolic Case

The first equation is the same as for the interface between the 2D-doubly parabolic and 1D-doubly parabolic case. Hence, we derive the same transmission condition reported in (42) for $u_{N_x+1,j}^{n+1}$ for $j=j_a,\dots ,j_b$.

For the flux, the transmission condition (17) gives us

$$v_0^{c,n+1}=-K(b-a)u_0^{c,n+1}+K▵y\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right),$$

(45)

This can be computed explicitly with the numerically computed values of $u_0^{c,n+1}$ and $u_{N_x+1,j}^{n+1}$.

Then, imposing the mass conservation:

$${\mathcal{I}}_{2\mathrm{D}}^{n+1}+{\mathcal{I}}_{2\mathrm{D}}^{n+1}-{\mathcal{I}}_{1\mathrm{D}}^n-{\mathcal{I}}_{1\mathrm{D}}^n=0$$

we get:

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{▵x}{2}\left[u_0^{c,n+1}-u_0^{c,n}+\lambda \frac{▵t}{▵x}\left(u_0^{c,n+1}-u_1^{c,n+1}\right)-\left(\frac{▵t}{▵x}-\frac{▵t}{2\lambda }\right)\left(-v_0^{c,n+1}-v_1^{c,n+1}\right)\right]\hfill \\ & +\frac{▵x▵t}{4\lambda }\left(F_{c,0}^{n+1}+F_{c,1}^{n+1}\right)\hfill \\ & +\frac{▵x▵y}{4}\left[{\displaystyle 4\sum _{j=j_a}^{j_b}\frac{▵tK}{▵x}\left(u_0^{c,n+1}+u_0^{c,n}-u_{N_x+1,j}^{n+1}-u_{N_x+1,j}^n\right)}\right]=0\hfill \end{array}\end{array}$$

We thus finally obtain the mass-preserving transmission condition, where we finally add the source term g, as:

$$\begin{array}{c}\hfill \begin{array}{ccc}{u}_0^{c,n+1}\hfill & =& u_0^{c,n}+\lambda \frac{▵t}{▵x}{\delta }_x^1{u}_0^{c,n+1}-\left(\frac{▵t}{▵x}-\frac{▵t}{2\lambda }\right)\left(v_0^{c,n+1}+v_1^{c,n+1}\right)\hfill \\ & & {\displaystyle -K\frac{▵t}{▵x}▵y\sum _{j=j_a}^{j_b}\left(u_0^{c,n+1}+u_0^{c,n}-u_{N_x+1,j}^{n+1}-u_{N_x+1,j}^n\right)}\hfill \\ & & -\frac{▵t}{2\lambda }\left(F_{c,0}^{n+1}+F_{c,1}^{n+1}\right)+\frac{▵t}{2}\left(g_0^{n+1}+g_1^{n+1}\right).\hfill \end{array}\end{array}$$

(46)

Proposition 2.

The complete numerical scheme derived for the 2D-doubly parabolic-1D-hyperbolic-parabolic model is mass-preserving—by construction—across the transmission conditions in absence of source terms.

Remark 5.

Note that the chemoattractant equation is the same as for the 1D-doubly parabolic and 2D-doubly parabolic case. Hence, the numerical schemes (33) and (19) with boundary conditions (32) and (43) can be used.

3.2.2. Multiple Channels

In the previous paragraphs, we connected the two-dimensional domain ${\Omega }_l$ with a single one-dimensional channel C at $(L_x,y)\in {\Omega }_l$ with $y\in \left[a_1,b_1\right]$, and $j_{a_1}$ and $j_{b_1}$, the positions of the endpoints of the corridor on the numerical grid. Of course, this can be extended to more channels.

Let $C_m$, with $m=1,\dots ,\mathcal{M}$ be $\mathcal{M}$ corridors, connected to the two-dimensional domain ${\Omega }_l$ at $\left(L_x,y_m\right)$ with $y_m\in \left[a_m,b_m\right]$ and $a_1>0$, $b_m<a_{m+1}$, for $m=1,\dots ,\mathcal{M}-1$, and $b_{\mathcal{M}}<L_y$ to avoid intersections of the corridors, with equal width $k▵y$, $k\in \mathbb{N}$.

3.3. Implemented Algorithm.

Before presenting the numerical tests in the next Section 4, we adapt 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 (CN) method on the diffusion and source term, which is a second-order implicit method and the explicit central method for the convection term. Moreover, since CN is only A-stable but not L-stable [38], we also need to choose a $\Delta t$ that is small enough to avoid spurious oscillations of the solution during transience.

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. Finally, the implicit-explicit numerical method used to compute the solutions for the density u in (8) inside the 2D domain ${\Omega }_l$ is:

$$\begin{array}{ccc}{u}_{i,j}^{n+1}\hfill & =& u_{i,j}^n+D_u\frac{▵t}{2}\left[\frac{{\delta }_x^2(u_{i,j}^n+u_{i,j}^{n+1})}{▵x^2}+\frac{{\delta }_y^2(u_{i,j}^n+u_{i,j}^{n+1})}{▵y^2}\right]\hfill \\ \hfill & & -\frac{▵t}{4}\left[\frac{{\delta }_x^0{F}_{i,j}^{x,n}}{▵x}+\frac{{\delta }_y^0{F}_{i,j}^{y,n}}{▵y}\right]+\frac{▵t}{2}\left(g_{i,j}^n+g_{i,j}^{n+1}\right)\hfill \\ \hfill & & \underset{\mathrm{artificial}\mathrm{viscosity}}{\underbrace{-▵t\left[\frac{{\delta }_x^2{\theta }_{i,j}^n}{2▵x}+\frac{{\delta }_y^2{\theta }_{i,j}^n}{2▵y}\right]}},\hfill \end{array}$$

(47)

with: ${\theta }_{i,j}^n:=\chi \left({\phi }_{i,k}^n\right)u_{i,j}^n|\nabla {\phi }_{i,j}^n|$ for $i=1,\dots ,N_x,j=1,\dots ,N_y$. As can be seen, the function $\theta$ used for the artificial viscosity is almost identical to f, with the exception of using the absolute value of $\nabla \phi$. 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, [36]. The numerical transmission condition on the left of node $1L$ ($i=N_x+1,j=j_a,\cdots ,j_b$) is:

$$\begin{array}{ccc}{u}_{N_x+1,j}^{n+1}\hfill & =& u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1(u_{N_x,j}^n+u_{N_x,j}^{n+1})\hfill \\ & & +D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)\hfill \\ & & +\frac{▵t}{▵x}\left(F_{N_x+1,j}^{x,n}+F_{N_x,j}^{x,n}\right)-\frac{▵t}{▵y}\left(F_{N_x+1,j+1}^{y,n}-F_{N_x+1,j-1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{N_x+1,j}^n+g_{N_x+1,j}^{n+1}\right)-▵t\left(\frac{{\delta }_x^1{\theta }_{N_x,j}^n}{▵x}+\frac{{\delta }_y^2{\theta }_{N_{x+1},j}^n}{2▵y}\right)\hfill \\ & & +\stackrel{(\mathrm{additional}\mathrm{term}\mathrm{for}\mathrm{transmission}\mathrm{condition})}{\underbrace{K\frac{▵t}{▵x}\left(u_0^{c,n}-u_{N_x+1,j}^n+u_0^{c,n+1}-u_{N_x+1,j}^{n+1}\right)}}.\hfill \end{array}$$

(48)

The role of permeability coefficient K in the positivity of (48) is discussed in Section 3.3.1.

For the external boundaries (the edges of the chamber ${\Omega }_l$ except at the junctions $j=j_a,\dots ,j_b$), we use:

$$\left\{\begin{array}{ccc}{u}_{i,0}^{n+1}\hfill & =\hfill & u_{i,0}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^2\left(u_{i,0}^n+u_{i,0}^{n+1}\right)+D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{i,0}^n+u_{i,0}^{n+1}\right)\hfill \\ & & -\frac{▵t}{2▵x}{\delta }_x^0{F}_{i,0}^{x,n}-\frac{▵t}{▵y}\left(F_{i,0}^{y,n}+F_{i,1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{i,0}^n+g_{i,0}^{n+1}\right)-▵t\left(\frac{{\delta }_x^2{\theta }_{i,0}^n}{2▵x}+\frac{{\delta }_y^1{\theta }_{i,0}^n}{▵y}\right),i=1,\dots ,N_x,\hfill \\ u_{i,N_y+1}^{n+1}\hfill & =\hfill & u_{i,N_y+1}^n+D_u\frac{▵t}{2▵x^2}{\delta }_x^2\left(u_{i,N_y+1}^n+u_{i,N_y+1}^{n+1}\right)\hfill \\ & & -D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{i,N_y}^n+u_{i,N_y}^{n+1}\right)\hfill \\ & & -\frac{▵t}{2▵x}{\delta }_x^0{F}_{i,N_y+1}^{x,n}+\frac{▵t}{▵y}\left(F_{i,N_y}^{y,n}+F_{i,N_y+1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{i,N_y+1}^n+g_{i,N_y+1}^{n+1}\right)-▵t\left(\frac{{\delta }_x^2{\theta }_{i,N_y+1}^n}{2▵x}+\frac{{\delta }_y^1{\theta }_{i,N_y}^n}{▵y}\right),i=1,\dots ,N_x,\hfill \\ u_{0,j}^{n+1}\hfill & =\hfill & u_{0,j}^n+D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{0,j}^n+u_{0,j}^{n+1}\right)+D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{0,j}^n+u_{0,j}^{n+1}\right)\hfill \\ & & -\frac{▵t}{▵x}\left(F_{0,j}^{x,n}+F_{1,j}^{x,n}\right)-\frac{▵t}{2▵y}{\delta }_y^0{F}_{0,j}^{y,n}\hfill \\ & & +\frac{▵t}{2}\left(g_{0,j}^n+g_{0,j}^{n+1}\right)-▵t\left(\frac{{\delta }_x^1{\theta }_{0,j}^n}{▵x}+\frac{{\delta }_y^2{\theta }_{0,j}^n}{2▵y}\right),j=1,\dots ,N_y,\hfill \\ u_{N_x+1,j}^{n+1}\hfill & =\hfill & u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,j}^n+u_{N_x,j}^{n+1}\right)\hfill \\ & & +D\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)\hfill \\ & & -\frac{▵t}{2▵x}{\delta }_x^0{F}_{i,N_y+1}^{x,n}+\frac{▵t}{▵x}\left(F_{i,N_y}^{y,n}+F_{i,N_y+1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{N_x+1,j}^n+g_{N_x+1,j}^{n+1}\right)-▵t\left(\frac{{\delta }_x^1{\theta }_{N_x,j}^n}{▵x}+\frac{{\delta }_y^2{\theta }_{N_x+1,j}^n}{2▵y}\right),j=1,\dots ,N_y,\hfill \\ & & j\ne j_a,\dots ,j_b.\hfill \end{array}\right.$$

(49)

For the corners, we use the following boundary conditions:

$$\left\{\begin{array}{ccc}{u}_{0,0}^{n+1}\hfill & =\hfill & u_{0,0}^n+D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{0,0}^n+u_{0,0}^{n+1}\right)+D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{0,0}^n+u_{0,0}^{n+1}\right)\hfill \\ & & -\frac{▵t}{▵x}\left(F_{0,0}^{x,n}+F_{1,0}^{x,n}\right)-\frac{▵t}{▵y}\left(F_{0,0}^{y,n}+F_{0,1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{0,0}^n+g_{0,0}^{n+1}\right)\hfill \\ & & -▵t\left(\frac{{\delta }_x^1{\theta }_{0,0}^n}{▵x}+\frac{{\delta }_y^0{\theta }_{0,0}^n}{▵y}\right),\hfill \\ u_{N_x+1,0}^{n+1}\hfill & =\hfill & u_{N_x+1,0}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,0}^n+u_{N_x,0}^{n+1}\right)\hfill \\ & & +D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{N_x+1,0}^n+u_{N_x+1,0}^{n+1}\right)\hfill \\ & & -\frac{▵t}{▵x}\left(F_{N_x+1,0}^{x,n}+F_{N_x,0}^{x,n}\right)-\frac{▵t}{▵y}\left(F_{N_x+1,0}^{y,n}+F_{N_x+1,1}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{N_x+1,0}^n+g_{N_x+1,0}^{n+1}\right)\hfill \\ & & -▵t\left(\frac{{\delta }_x^1{\theta }_{N_x,0}^n}{▵x}+\frac{{\delta }_y^1{\theta }_{N_x+1,0}^n}{▵y}\right),\hfill \\ u_{0,N_y+1}^{n+1}\hfill & =\hfill & u_{0,N_y+1}^n+D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{0,N_y+1}^n+u_{0,N_y+1}^{n+1}\right)\hfill \\ & & -D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{0,N_y}^n+u_{0,N_y}^{n+1}\right)\hfill \\ & & -\frac{▵t}{▵x}\left(F_{0,N_y+1}^{x,n}+F_{1,N_y+1}^{x,n}\right)-\frac{▵t}{▵y}\left(F_{0,N_y+1}^{y,n}+F_{0,N_y}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{0,N_y+1}^n+g_{0,N_y+1}^{n+1}\right)\hfill \\ & & -▵t\left(\frac{{\delta }_x^1{\theta }_{0,N_y+1}^n}{▵x}+\frac{{\delta }_y^1{\theta }_{0,N_y}^n}{▵y}\right),\hfill \\ u_{N_x+1,N_y+1}^{n+1}\hfill & =\hfill & u_{N_x+1,N_y+1}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1\left(u_{N_x,N_y+1}^n+u_{N_x,N_y+1}^{n+1}\right)\hfill \\ & & -D_u\frac{▵t}{▵y^2}{\delta }_y^1\left(u_{N_x+1,N_y}^n+u_{N_x+1,N_y}^{n+1}\right)\hfill \\ & & -\frac{▵t}{▵x}\left(F_{N_x+1,N_y+1}^{x,n}+F_{N_x,N_y+1}^{x,n}\right)\hfill \\ & & -\frac{▵t}{▵y}\left(F_{N_x+1,N_y+1}^{y,n}+F_{N_x+1,N_y}^{y,n}\right)\hfill \\ & & +\frac{▵t}{2}\left(g_{N_x+1,N_y+1}^n+g_{N_x+1,N_y+1}^{n+1}\right)\hfill \\ & & -▵t\left(\frac{{\delta }_x^1{\theta }_{N_x,N_y+1}^n}{▵x}+\frac{{\delta }_y^1{\theta }_{N_x+1,N_y}^n}{▵y}\right).\hfill \end{array}\right.$$

(50)

Similarly, for the chemoattractant $\varphi$, we have the implicit-explicit scheme in the interior points of the 2D domain:

$$\begin{array}{ccc}{\varphi }_{i,j}^{n+1}\hfill & =& {\varphi }_{i,j}^n+D_{\varphi }\frac{▵t}{2}\left[\frac{{\delta }_x^2\left({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1}\right)}{▵x^2}+\frac{{\delta }_y^2\left({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1}\right)}{▵y^2}\right]\hfill \\ \hfill & & \frac{▵t}{2}(a(u_{i,j}^n+u_{i,j}^{n+1})-\frac{▵t}{2}(b({\varphi }_{i,j}^n+{\varphi }_{i,j}^{n+1})\hfill \end{array}$$

(51)

At the boundaries and the corners of the numerical schemes for $\varphi$, we use, respectively, conditions:

$$\left\{\begin{array}{ccc}{\varphi }_{0,0}^{n+1}\hfill & =\hfill & {\varphi }_{0,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{0,0}^n+{\varphi }_{0,0}^{n+1})+2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{0,0}^n+{\varphi }_{0,0}^{n+1})\hfill \\ \hfill & \hfill & +a▵t(u_{0,0}^n+u_{0,0}^{n+1})-b▵t({\varphi }_{0,0}^n+{\varphi }_{0,0}^{n+1}),\hfill \\ {\varphi }_{N_x+1,0}^{n+1}\hfill & =\hfill & {\varphi }_{N_x+1,0}^n-D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x+1,0}^n+{\varphi }_{N_x+1,0}^{n+1})\hfill \\ \hfill & \hfill & +D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{N_x+1,0}^n+{\varphi }_{N_x+1,0}^n)\hfill \\ \hfill & \hfill & +a▵t(u_{N_x+1,0}^n+u_{N_x+1,0}^{n+1})-b▵t({\varphi }_{N_x+1,0}^n+{\varphi }_{N_x+1,0}^{n+1}),\hfill \\ {\varphi }_{N_x+1,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{Nx+1,N_y+1}^n-D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x,N_y+1}^n+{\varphi }_{N_x,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -D_{\varphi }\frac{▵t}{▵y^2}{\delta }_x^1({\varphi }_{N_x+1,N_y+1}^n+{\varphi }_{N_x+1,N_y}^n)\hfill \\ \hfill & \hfill & +a▵t{u}_{N_x+1,N_y+1}^n-b▵t{\varphi }_{N_x+1,N_y+1}^n,\hfill \\ {\varphi }_{0,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{0,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{1,N_y+1}^n+{\varphi }_{1,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{0,N_y}^n+{\varphi }_{0,N_y}^{n+1})+a▵t{u}_{0,N_y+1}^n-b▵t{\varphi }_{0,N_y+1}^n,\hfill \end{array}\right.$$

(52)

For the borders, we use:

$$\left\{\begin{array}{ccc}{\varphi }_{i,0}^{n+1}\hfill & =\hfill & {\varphi }_{i,0}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^2({\varphi }_{i,0}^n+{\varphi }_{i,0}^{n+1})+2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{i,0}^n+{\varphi }_{i,0}^{n+1})\hfill \\ \hfill & \hfill & +a▵t{u}_{i,0}^n-b▵t{\varphi }_{i,0}^n,\hfill \\ {\varphi }_{i,N_y+1}^{n+1}\hfill & =\hfill & {\varphi }_{i,N_y+1}^n+D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^2({\varphi }_{i,N_y+1}^n+{\varphi }_{i,N_y+1}^{n+1})\hfill \\ \hfill & \hfill & -2D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^1({\varphi }_{i,N_y}^n+{\varphi }_{i,N_y}^{n+1})+a▵t{u}_{i,N_y+1}^n-b▵t{\varphi }_{i,N_y+1}^n,\hfill \\ {\varphi }_{0,j}^{n+1}\hfill & =\hfill & {\varphi }_{0,j}^n+D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^2({\varphi }_{0,j}^n+{\varphi }_{0,j}^{n+1})+2D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{0,j}^n+{\varphi }_{0,j}^{n+1})\hfill \\ \hfill & \hfill & +a▵t{u}_{0,j}^n-b▵t{\varphi }_{0,j}^n,\hfill \\ {\varphi }_{N_x+1,j}^{n+1}\hfill & =\hfill & {\varphi }_{N_x+1,j}^n+D_{\varphi }\frac{▵t}{▵y^2}{\delta }_y^2({\varphi }_{N_x+1,j}^n+{\varphi }_{N_x+1,j}^{n+1})\hfill \\ \hfill & \hfill & -2D_{\varphi }\frac{▵t}{▵x^2}{\delta }_x^1({\varphi }_{N_x,j}^n+{\varphi }_{N_x,j}^{n+1})+a▵t{u}_{N_x+1,j}^n-b▵t{\varphi }_{N_x+1,j}^n.\hfill \end{array}\right.$$

(53)

Note that for $i=N_x+1$ the last formula in (53) is applied for $j=1,\dots ,N_y,j\ne j_a,\dots ,j_b$.

Remark 6.

If we consider a two-dimensional domain ${\Omega }_r$ connected to the right endpoint of the one-dimensional corridor C, the complete numerical scheme for the left domain ${\Omega }_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 ${\Omega }_l$ and the right for the corridor) are reversed to the left for the corridor and the right for the box ${\Omega }_r$. In the numerical scheme, the change only affects the channel C, where we have transmission conditions also for $u_{N+1}^{c,n}$ (resp. $v_{N+1}^{c,n}$). 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_0^{c,n}$ (resp. $v_0^{c,n}$).

For the computation of solutions on the one-dimensional channel C, we have two different approximations depending on the choice of the model we assign to it. If we solve the doubly parabolic problem (9), the approximation scheme used is the Crank–Nicolson scheme, as above:

$$\begin{array}{ccc}{u}_i^{c,n+1}\hfill & =& u_i^n+D_{u_c}\frac{▵t}{2}\left[\frac{{\delta }_x^2\left(u_i^{c,n}+u_i^{c,n+1}\right)}{▵x^2}\right]-▵t{\delta }_x^0{F}_{c,i}^n\hfill \\ \hfill & & +\frac{▵t}{2}\left(g_i^n+g_i^{n+1}\right)-▵t\left(\frac{{\delta }_x^2{\theta }_i^n}{2▵x}\right),\hfill \end{array}$$

(54)

with the transmission condition on the left of node $1L$ ($i=0$) given by:

$$\begin{array}{ccc}{u}_0^{c,n+1}\hfill & =& u_0^{c,n}+D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_0^{c,n}+u_0^{c,n+1}\right)-\frac{▵t}{▵x}\left(F_{c,0}^n+F_{c,1}^n\right)-▵t\frac{{\delta }_x^1{\theta }_0^n}{▵x}\hfill \\ \hfill & & -K\frac{▵t}{▵x}\left\{{\displaystyle (b-a)(u_0^{c,n}+u_0^{c,n+1})+\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)}\right\}\hfill \\ \hfill & & +\frac{▵t}{2}\left(g_0^n+g_0^{n+1}\right).\hfill \end{array}$$

(55)

If, instead, we need to solve the hyperbolic-parabolic problem (7), an implicit version of the second-order AHO scheme is applied—see the scheme reported in (44)—in order to ensure the stability of numerical solutions in the channels.

The scheme (44) is endowed with transmission conditions (45) and (46).

3.3.1. Stability at Interfaces

Note that, in order to ensure the positivity of the quantities in the above formulas deriving from the KK conditions—i.e., (48) for the 2D domain and (55) or (45) 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 (48) and (45), one needs to ensure that $K\frac{▵t}{▵x}$ and, respectively, $K\frac{▵t}{▵x}▵y$ is not too big in order to damp possible high oscillations produced by the term in parentheses. Similarly, in (55) we need to check that $K\frac{▵t}{▵x}(b-a)$ is small in order to prevent the growth of the negative term.

Moreover, as previously discussed, we need to check that the numerical monotonicity condition is satisfied:

$$\begin{array}{c}\hfill \frac{k_1}{{\left(k_2+{\phi }_{i,j}^n\right)}^2}\left|{\partial }_{x,i,j}^n{\phi }_{i,j}^n\right|\le \sqrt{D_M}\end{array}$$

(56)

in the computational domain in order to ensure non-negative solutions.

Now, we consider the interface between the 2D and 1D domains. If we assume $g=0$ and $f=\chi \left(\varphi \right)\nabla \varphi$, the first equation in the 2D parabolic system, (8), rewrites as:

$${\partial }_{t}u=D▵u-div\left(u·f\right).$$

(57)

The transmission condition (48) reads as:

$$\begin{array}{ccc}{u}_{N_x+1,j}^{n+1}\hfill & =\hfill & u_{N_x+1,j}^n-D_u\frac{▵t}{▵x^2}{\delta }_x^1(u_{N_x,j}^n+u_{N_x,j}^{n+1})+D_u\frac{▵t}{2▵y^2}{\delta }_y^2\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)\hfill \\ & & +\frac{▵t}{▵x}\left(f_{N_x,j}^{x,n}{u}_{N_x,j}^n+f_{N_x+1,j}^{x,n}{u}_{N_x+1,j}^n\right)\hfill \\ & & -\frac{▵t}{2▵y}\left(f_{N_x+1,j+1}^{y,n}{u}_{N_x+1,j+1}^n-f_{N_x+1,j-1}^{y,n}{u}_{N_x+1,j-1}\right)\hfill \\ & & +\frac{▵t}{▵x}\left(|f_{N_x,j}^{x,n}|u_{N_x,j}^n-\left|f_{N_x+1,j}^{x,n}\right|u_{N_x+1,j}^n\right)\hfill \\ & & +\frac{▵t}{2▵y}\left(|f_{N_x+1,j+1}^{y,n}|u_{N_x+1,j+1}^n-2\left|f_{N_x+1,j}^{y,n}\right|u_{N_x+1,j}^n\right.\hfill \\ & & +\left.|f_{N_x+1,j-1}^{y,n}|u_{N_x+1,j-1}^n\right)\hfill \\ & & +\underset{\mathrm{KK}-\mathrm{transmission}\mathrm{term}}{\underbrace{\frac{▵x}{▵t}K\left(\left(u_0^{c,n}-v_{N_x+1,j}^n\right)+\left(u_0^{c,n+1}-u_{N_x+1,j}^{n+1}\right)\right)}}.\hfill \end{array}$$

(58)

Then, the transmission condition (55) for $g=0$ and $f=\chi \left(\varphi \right){\varphi }_x$ now reads as:

$$\begin{array}{ccc}{u}_0^{c,n+1}\hfill & =\hfill & u_0^{c,n}+D_{u_c}\frac{▵t}{▵x^2}{\delta }_x^1\left(u_0^{c,n}+u_0^{c,n+1}\right)-\frac{▵t}{▵x}\left(u_0^{c,n}{f}_0^n+u_1^{c,n}{f}_1^n\right)\hfill \\ & & +\frac{▵t}{▵x}\left(u_1^{c,n}|f_1^n|-u_0^{c,n}\left|f_0^n\right|\right)\hfill \\ & & +\underset{\mathrm{KK}-\mathrm{transmission}\mathrm{term}}{\underbrace{\frac{▵x}{▵t}K\left({\displaystyle \sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)-(b-a)\left(u_0^n+u_0^{n+1}\right)}\right)}}\hfill \end{array}$$

(59)

for $j=j_a,\dots j_b$. For the transmission condition (58), we see that monotonicity is preserved when:

$$\begin{array}{c}1-D_u\frac{▵t}{▵x^2}-D_u\frac{▵t}{▵y^2}+\frac{▵t}{▵x}{f}_{N_x+1,j}^{x,n}-\frac{▵t}{▵x}|f_{N_x+1,j}^{x,n}|-\frac{▵t}{▵y}\left|f_{N_y+1,j}^{y,n}\right|-\frac{▵t}{▵x}K>0,\hfill \end{array}$$

(60)

which gives us the stability condition for the left side of the interface.

For the right side of the interface, if we assign the doubly parabolic 1D system (9), we have the stability condition:

$$\begin{array}{cc}& 1-D_{u_c}\frac{▵t}{▵x^2}+\lambda f_{c,0}^n-\frac{▵t}{▵x}\left|f_0^n\right|-\frac{▵t}{▵x}K(b-a).\hfill \end{array}$$

(61)

We underline that (61) is not only influenced by the KK-constant K but also by the channel width $\sigma :=(b-a)$, which must be taken care of accordingly.

Analogously, we conduct the derivation of stability conditions for the transmission conditions in the case where we assign the 2D parabolic model (57) on the left part of the interface and the hyperbolic system (10) on the right. For the sake of clarity, we rewrite the hyperbolic part of the system (10) for the density flux of individuals when we assume $g=0$:

$$\left\{\begin{array}{ccc}{\partial }_t{u}_c+{\partial }_x{v}_c\hfill & =& 0,\hfill \\ {\partial }_t{v}_c+{\lambda }_c^2{\partial }_x{u}_c\hfill & =& F_c.\hfill \end{array}\right.$$

(62)

The KK-transmission explicit version of conditions (46) and (45) in this case read as:

$$\begin{array}{c}{u}_0^{c,n+1}=u_0^{c,n}+\lambda \frac{▵t}{▵x}\left(u_1^{c,n}-u_0^{c,n}\right)-\left(\frac{▵t}{▵x}-\frac{▵t}{2\lambda }\right)\left(v_0^{c,n}+v_1^{c,n}\right)-\frac{▵t}{2\lambda }\left(F_{c,0}^n+F_{c,1}^n\right)\hfill \\ {\displaystyle -K\frac{▵t}{▵x}▵y\sum _{j=j_a}^{j_b}\left(u_0^{c,n}-u_{N_x+1,j}^n+u_0^{c,n+1}-u_{N_x+1,j}^{n+1}\right),}\hfill \\ u_{N+1}^{c,n+1}=u_{N+1}^{c,n}+\lambda \frac{▵t}{▵x}\left(u_N^{c,n}-u_{N+1}^{c,n}\right)+\left(\frac{▵t}{▵x}-\frac{▵t}{2\lambda }\right)\left(v_N^{c,n}+v_{N+1}^{c,n}\right)+\frac{▵t}{2\lambda }\left(F_{c,N}^n+F_{c,N+1}^n\right)\hfill \\ {\displaystyle -K\frac{▵t}{▵x}▵y\sum _{j=j_a}^{j_b}\left(u_{N+1}^{c,n}-u_{0,j}^{c,n}+u_{N+1}^{c,n+1}-u_{0,j}^{c,n+1}\right),}\hfill \\ v_0^{c,n}=-K(b-a)u_0^{c,n}+K▵y\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^{l,n}+u_{N_x+1,j}^{l,n+1}\right),\hfill \\ v_{N+1}^{c,n}=-K(b-a)u_{N+1}^{c,n}+K▵y\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^{r,n}+u_{N_x+1,j}^{r,n+1}\right).\hfill \end{array}$$

(63)

In order to establish the monotonicity condition, as in [39], we need to diagonalize the boundary conditions above into the diagonal variables $w_0^{n+1},z_0^{n+1},w_{N+1}^{n+1},z_{N+1}^{n+1}$ with the relation $u=w+z$ and $v=\lambda \left(z-w\right)$.

We have:

$$\begin{array}{c}\hfill \begin{array}{ccc}{v}_0^{c,n+1}\hfill & =& \lambda \left(z_0^{c,n+1}-w_0^{c,n+1}\right)=-K(b-a)(z_0^{c,n+1}+w_0^{c,n+1})\hfill \\ & & +K▵y\sum _{j=j_a}^{j_b}\left(u_{N_x+1,j}^{l,n+2}+u_{N_x+1,j}^{l,n+1}\right).\hfill \end{array}\end{array}$$

If we set $\rho :=\frac{\lambda -(b-a)K}{\lambda +(b-a)K}$,${\varsigma }^{n+1}:=+\frac{K}{\lambda +(b-a)K}{\sum }_{j=j_a}^{j_b}\left(u_{N_x+1,j}^{l,n+2}+u_{N_x+1,j}^{l,n+1}\right)$ we find

$$z_0^{n+1}=\rho w_0^{n+1}+▵y{\varsigma }^{n+1}$$

and

$$\begin{array}{c}\hfill \begin{array}{ccc}{w}_0^{n+1}+\rho w_0^{n+1}+{\varsigma }^{n+1}\hfill & =& (1+\rho )w_0^n+{\varsigma }^n+2\lambda \frac{▵t}{▵x}\left(w_1^n-\rho w_0^n-{\varsigma }^n\right)\hfill \\ & & +\frac{▵t}{2}\left(\left(\rho -1\right)w_0^n+{\varsigma }^n+z_1^n-w_1^n\right)-\frac{▵t}{2\lambda }\left(F_{c,0}^n+F_{c,1}^n\right)\hfill \\ & & -\frac{▵t}{▵x}{K}_k((b-a)\left(1+\rho \right)\left(w_0^n+w_0^{n+1}\right)+(b-a)\left({\varsigma }^n+{\varsigma }^{n+1}\right)\hfill \\ & & -\sum _{j=j_{a_k}}^{j_{b_k}}\left(u_{N_x+1,j}^n+u_{N_x+1,j}^{n+1}\right)).\hfill \end{array}\end{array}$$

Now, by applying the monotonicity condition we get the following inequality:

$$\begin{array}{c}\left(1+\rho \right)-2\lambda \frac{▵t}{▵x}\rho +\frac{▵t}{2}\left(\rho -1\right)-\frac{▵t}{▵x}K(b-a)\left(1+\rho \right)-\frac{▵t}{2\lambda }{F}_{c,0}^n>0\hfill \\ \iff ▵t\le \frac{1+\rho }{\frac{2\lambda }{▵x}\rho -\frac{\rho -1}{2}+K\frac{(b-a)(1+\rho )}{▵x}+\frac{F_{c,0}^n}{2\lambda }}\hfill \end{array}$$

(64)

For the implicit AHO, the condition above reads as:

$$▵t\le \frac{1+\rho }{K\frac{(b-a)(1+\rho )}{▵x}+\frac{F_{c,0}^n}{2\lambda }}.$$

In Figure 3, the time step restriction (64) is depicted for a qualitative understanding of the effect of the Kedem–Katchalsky constant K and of the channel width $\sigma$ on the time step $▵t$. As expected, the time step $▵t$ must be chosen smaller when either K or $\sigma$ increases. Furthermore for $K=0$ we recover the time step restriction of the $AH{O}^2$-scheme $▵t\le \frac{4▵x}{▵x+4\lambda }=2·{10}^{-3}$. Since the values of K are typically of similar magnitude to the diffusion coefficients, the additional stability restrictions caused by the hyperbolic part of the transmission conditions are minimal.

Finally, we also point out that the time step restriction for the transmission condition for the one-dimensional parabolic Equation (61) is much more severe than for the one-dimensional hyperbolic (64), which can also be seen qualitatively in the steepness of Figure 3 and Figure 4.

Comparing all time step restrictions (61) and (64) with each other, it is evident that the restriction for the two-dimensional parabolic transmission condition dominates the full model.

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 [40] can be used, which take advantage of the mostly sparse structure of the Jacobian matrix.

4. Numerical Tests and Results

We start this section with a preliminary test on mass-preserving properties at the boundaries. Indeed, in absence of source terms, the masses of cells and chemical substances are preserved. Then, in order to perform a numerical verification of this property, we consider the numerical approximation at the interface between 1D-1D models in the next numerical example.

Example 1.

In Figure 5, a comparison between the central mass-preserving (24) and standard finite difference (23) boundary conditions is depicted, 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 (on the right). From this 1D numerical example, the necessity of developing modified boundary conditions which are consistent and preserve the mass correctly is evident. We also underline that for the discretization of the Neumann condition, the forward scheme is first-order accurate, while the central scheme is second-order accurate; for more details, see [41].

From now on, this section is devoted to the presentation of the numerical tests and the parameters of the problem are reported in

Table 1. Parameters of the problem.

Parameter	Description	Units	Value	Ref.
$D_M$	Diffusivity of cells	$\mathsf{\mu }$m${}^2$/s	$9\times {10}^2$	[33]
$D_T$	Diffusivity of T -scenario 2	$\mathsf{\mu }$m${}^2$/s	$56\times {10}^1$	empirical
$D_T$	Diffusivity of T -scenario 1	$\mathsf{\mu }$m${}^2$/s	$5.6\times {10}^1$	[33]
$D_{\phi },D_{\omega }$	Diffusivity of chemoattractants	$\mathsf{\mu }$m${}^2$/s	$2\times {10}^2$	[33]
${\alpha }_{\phi }$	growth rate of $\phi$ in scenario 2	s${}^{-1}$/cell	0	empirical
${\alpha }_{\phi }$	growth rate of $\phi$ for scenario 1	s${}^{-1}$/cell	$0.5\times {10}^{-1}$	[42]
${\beta }_{\phi }$	consumption rate of $\phi$	s${}^{-1}$	${10}^{-4}$	[42]
${\alpha }_{\omega }$	growth rate of $\omega$ for scenario 1	s${}^{-1}$/cell	${10}^{-1}$	[42]
${\alpha }_{\omega }$	growth rate of $\omega$ for scenario 2	s${}^{-1}$/cell	0	empirical
${\beta }_{\omega }$	consumption rate of $\omega$	s${}^{-1}$	${10}^{-4}$	[42]
$k_1$	cellular drift velocity	mol cm${}^2$s${}^{-1}$	$3.9\times {10}^{-9}$	[33]
$k_2$	receptor dissociation constant	mol	$5\times {10}^{-6}$	[33]
S	maximum secretion rate of the chemicals	g/$\mathsf{\mu }$m${}^3$	$1.9676\times {10}^{-11}$	[33]
$\gamma$	equivalent Michaelis constant	cells$/\mathsf{\mu }$m${}^3$	${10}^{-4}$	[33]
L	length of the corridor	$\mathsf{\mu }$m	500	datum
$L_x$	horizontal size of the box	$\mathsf{\mu }$m	100	datum
$L_y$	vertical size of the box	$\mathsf{\mu }$m	1000	datum
K	Kedem–Katchalsky constant	-	500	empirical

. 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, [1,2,25].

We remark that we decided to perform numerical simulations of the chip geometry by assigning the 1D-hyperbolic-parabolic model on channels, since it seems more realistic.

Example 2.

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 ${\Omega }_l$ with one corridor positioned in the middle and only one cell family with initial distribution $u(x,y,0)=5e^{-\frac{1}{2}\left({(x-0.5)}^2+{(y-0.5)}^2\right)}$. Since we do not have any analytical solution for this problem, we choose $\Delta t$ and $\Delta x=\Delta y$, which are small enough to obtain reasonable error estimations. In this case, we use $dt={10}^{-4}$ and $\Delta x=\Delta y=5\times {10}^{-4}$ for the approximation $u_e$ at time $t=100$ and calculate the error as the quantity $\parallel u_e-u_{\mathrm{approx}}\parallel$ in $L^1$-norm.

In order to confirm the order of our scheme, we use a log-log-plot with constant and small enough $\Delta t$ (resp. $\Delta x$) and decreasing $\Delta x$ (resp. $\Delta t$). As shown in Figure 6 and Figure 7, the time order and space order are equal to line with slope 2 in the log-log plot, which corresponds to our scheme of order 2 in time and space.

Now, we describe the simulation of the chip environment. All the simulations were performed in MATLAB©. The computational time for a simulation on the complete geometry until time t = 50,000 took about 400 s on an Intel(R) Core(TM) i7-3630 QM CPU 2.4 GHz. The computational domain is schematized in Figure 2 with the two chambers and 5 corridors $C_m:=[0,L],m=1,\dots ,5$ with the same width $=12\mathsf{\mu }$m and equispaced from each other. The numerical method implemented is listed in Section 3.3; for the 1D channels, the $AH{O}^2$-scheme (of second order) is implemented, since there we are considering the hyperbolic-parabolic model. The discretization grid has time step size $▵t=100$ s and space size $▵x=2.5\mathsf{\mu }\mathrm{m},▵y=25\mathsf{\mu }$m.

For the examples below, we assume the following initial condition (time $t=0$) for the tumor cells distribution on the chip for $(x,y)\in {\Omega }_l$:

$$T(x,y,0)=5e^{-{10}^{-4}\left(x^2+{(y-500)}^2\right)}+5e^{-{10}^{-4}\left(x^2+{(y-5)}^2\right)}+5e^{-{10}^{-4}\left(x^2+{(y-1000)}^2\right)},$$

(65)

Whereas, in the corridors and the right chamber, no tumor cells are present.

For the immune cells distribution on the chip for $(x,y)\in {\Omega }_r$, we assign:

$$M(x,y,0)=5,\mathrm{for}x,y\in {\Omega }_r,$$

(66)

to mean that macrophages are disposed in the right chamber, whereas no immune cells are present in the left chamber nor in the corridors at the beginning.

For the chemoattractants, we set a constantly null initial density for $\omega$ and $\phi$ in both the chambers and also in the channels.

Example 3.

First scenario: treated case. For the following numerical simulation, we replicate the laboratory experiment of the first scenario (treated case) described in Section 2.1; see also [1]. All the parameters are reported in Table 1.

The results are depicted in the following Figure 8 and Figure 9.

In Figure 8 and Figure 9, we can see the density of the tumor cells T and macrophages M for different times $t=0,t=10000 s$, and t = 50,000 s. 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 tumor cells are previously treated by a chemoterapy drug, they slowly diffuse around and stay confined in the left chamber ${\Omega }_l$ during all the simulation time; in the meanwhile they produce chemoattractant φ attracting immune cells. Immune cells M, instead, diffuse around in ${\Omega }_r$, cross the channels, and after a certain time they enter the left chamber ${\Omega }_l$ while creating chemoattractant ω.

This is due to the fact that the chemoattractant φ produced by cancer cells travels through channels and induces 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, as we observe from the laboratory experiment.

At the final time t = 50,000 s, we can see that the quantity of tumor cells is decreasing under the action of immune cells producing chemokine ω.

Example 4.

Second scenario: untreated case.In this numerical test, we consider the second scenario where the tumor is not treated with any medicine. Therefore, in this case we assume a higher diffusion coefficient for the tumor, but the initial conditions are the same as those used above. The results are depicted in the following Figure 10.

In Figure 10, we can see the density of the tumor cells T and macrophages M for times t = 10,000 s and t = 50,000 s. 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, in the laboratory experiment, untreated tumor cells diffuse around, cross the channels, and enter the right chamber ${\Omega }_r$ after some time, we try to reproduce such behavior by using a diffusion coefficient $D_T$ of a order of magnitude higher respect to the one used in the first scenario. Moreover, in this case the tumor cells do not produce chemoattractant φ; for this reason, here we set ${\alpha }_{\varphi }=0$. Immune cells M diffuse around in ${\Omega }_r$ and do not cross the channels, since the chemical stimulus is not secreted by T cells. Moreover, the production of the chemical substance w is neglected in this case, thus tumor cells are not killed.

We only mention that we tested the 1D-doubly parabolic model on channels and compared it with the hyperbolic-parabolic model used in the previous examples. By using the same initial data as for the other examples, we notice that the doubly parabolic model seems to have a similar pattern as for the hyperbolic-parabolic model depicted above, but the scale of the quantities differs a lot between these models.

In particular, for the doubly parabolic model the concentration of the tumor cells T is two or three order of magnitude higher. This is due to the much slower movement of the immune cells through the channels. This also explains the much higher concentration of the chemoattractant $\varphi$ because of the much higher concentration of T compared to in the other models.

In the following Figure 11, we represent the density of tumor cells and immune cells as particles, 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.

5. Conclusions 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 a 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 the 2D domain and the mass-preserving transmission conditions at the 2D–1D interfaces.

Furthermore, from the modelling point of view, we studied the dynamics in the channels in the case of a doubly parabolic model and a hyperbolic-parabolic model. Since we obtained comparable asymptotic states, we decided to apply the hyperbolic-parabolic model in order to obtain a finite speed of propagation in the channels, which seems to be more realistic. In this framework, bearing in mind the laboratory experiments on a chip described in Section 2.1, 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 obtained from cell tracking in a microfluidic chip [1,2,25].