Hybrid models of chemotaxis with application to leukocyte migration

Abstract

Many chemical and biological systems involve reacting species with vastly different numbers of molecules/agents. Hybrid simulations model such phenomena by combining discrete (e.g., agent-based) and continuous (e.g., partial differential equation- or PDE-based) descriptors of the dynamics of reactants with small and large numbers of molecules/agents, respectively. We present a stochastic hybrid algorithm to model a stage of the immune response to inflammation, during which leukocytes reach a pathogen via chemotaxis. While large numbers of chemoattractant molecules justify the use of a PDE-based model to describe the spatiotemporal evolution of its concentration, relatively small numbers of leukocytes and bacteria involved in the process undermine the veracity of their continuum treatment by masking the effects of stochasticity and have to be treated discretely. Motility and interactions between leukocytes and bacteria are modeled via random walk and a stochastic simulation algorithm, respectively. Since the latter assumes the reacting species to be well mixed, the discrete component of our hybrid algorithm deploys stochastic operator splitting, in which the sequence of the diffusion and reaction operations is determined autonomously during each simulation step. We conduct a series of numerical experiments to ascertain the accuracy and computational efficiency of our hybrid simulations and, then, to demonstrate the importance of randomness for predicting leukocyte migration and fate during the immune response to inflammation.

Appendices

GMP simulation

In each mesh element (i, j) at any time t, the system is described by the numbers of leukocytes and bacteria ${\mathbf {n}}^{(i,j)}(t)$, as defined in (12) and (13). The superscript (i, j) is omitted below to simplify the notation.

Let ${\mathbb {P}}_0[\tau | {\mathbf {n}} ,t]$ denote the (conditional) probability of no reactions taking place during the time interval $[t,t+\tau )$ provided that the system is at state ${\mathbf {n}}$ at time t. The reaction system is assumed to be Markovian, i.e., the probability that no reactions occur during $[t,t+\tau + \text{ d }\tau ]$ equals the product of the probabilities of no reactions occurring during $[t,t+\tau )$ and during $[t+\tau ,t+\tau + \text{ d }\tau )$. A propensity function $a_r$ is defined such that $a_r \text{ d }\tau$ is the probability that both the next reaction will be the rth reaction and it will occur during $[t+\tau ,t+\tau + \text{ d }\tau ]$. Then, one obtains (Gillespie 1976)

$$\begin{aligned} \begin{aligned}&{\mathbb {P}}_0[\tau + \text{ d }\tau | {\mathbf {n}}, t] = {\mathbb {P}}_0[\tau | {\mathbf {n}},t](1-a_{\text {sum}}({\mathbf {n}}) \text{ d }\tau ),\\&a_{\text {sum}}({\mathbf {n}}) = \sum _{r=1}^S a_r({\mathbf {n}}), \end{aligned} \end{aligned}$$

(25)

where S is the number of chemical reactions. Taking the limit as $\text{ d }\tau \rightarrow 0$ and solving the resulting PDE leads to

$$\begin{aligned} {\mathbb {P}}_0(\tau | {\mathbf {n}},t) = {\text {e}}^{-a_{\text {sum}}({\mathbf {n}})\tau }. \end{aligned}$$

(26)

The propensity functions for the reactions in (3) are listed in Table 3.

Table 3 Propensity functions for the leukocytes-bacteria-chemoattractant system

It follows from the definition of ${\mathbb {P}}_0$ and $a_r$ that the probability ${\mathbb {P}}(\tau ,r | {\mathbf {n}} ,t)$ of both the next reaction being the rth reaction and occurring during $[t+\tau ,t+\tau +\text{ d }\tau )$, given the present state of the system ${\mathbf {n}}(t)$, is

$$\begin{aligned} {\mathbb {P}}(\tau ,r| {\mathbf {n}} ,t) = {\mathbb {P}}_0[\tau | {\mathbf {n}},t]a_r({\mathbf {n}}). \end{aligned}$$

(27)

Accounting for (26),

$$\begin{aligned} {\mathbb {P}}(\tau ,r| {\mathbf {n}} ,t) = \frac{a_r( {\mathbf {n}})}{a_{\text {sum}}({\mathbf {n}})}a_{\text {sum}}({\mathbf {n}}) {\text {e}}^{-a_{\text {sum}}({\mathbf {n}})\tau }. \end{aligned}$$

(28)

The ratio $a_r({\mathbf {n}}) / a_{\text {sum}}({\mathbf {n}})$ represents the probability of a discrete random variable: the label of the next reaction. The term $a_{\text {sum}}({\mathbf {n}}) \exp [-a_{\text {sum}}({\mathbf {n}})\tau ]$ is the exponential density function of a continuous random variable: the time at which the next reaction will occur.

The GMP algorithm for reactions is presented below.

Additional numerical tests

1.1 Particle-based solution of a PDE

Consider a one-dimensional advection-diffusion PDE on the semi-infinite domain ${\varOmega }= [0,+\infty )$,

$$\begin{aligned} \left\{ \begin{aligned}&\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}=D\frac{\partial ^2 u}{\partial x^2},\quad u(x,t_0) = u_0(x),\\&vu(0,t)+D\frac{\partial u}{\partial x}(0,t) =0,\quad u(+\infty ,t)=0, \end{aligned} \right. \end{aligned}$$

(29)

where $u_0$ is non-zero on a single element of the numerical mesh. This problem can be represented by a stochastic equation with a reflective boundary condition at $x = 0$:

$$\begin{aligned} X^k(t+{\varDelta }t) = X^k(t)+ v{\varDelta }t + \sqrt{2D{\varDelta }t}\xi , \end{aligned}$$

(30)

where $k=1,\cdots ,M$ and $\xi$ is a standard Gaussian random variable. $X^k$ is the position of particle k, driven by Brownian dynamics and velocity v. Initially at $t_0 = 0$, the M particles are located very close to $x=0$ to mimic the Dirac delta point source, i.e.,

$$\begin{aligned} X^k(0) = \epsilon ,\quad {\text {where }}0 < \epsilon \ll 1. \end{aligned}$$

(31)

Fig. 8

Left column: Deterministic and stochastic solutions of (29). Results are shown at $t=500~\hbox {s}$, 1000 s, and 1500 s with Péclet numbers ${\text {Pe}}=0.5$, 2.0, and 10.0; Right column: Model discrepancy $\epsilon$ at $t=500~\hbox {s}$, 1000 s, and 1500 s. The error $\epsilon$ is plotted in log–log scale as function of the number of particles $M = 2^p\times 1000$ ($p=1,\ldots ,8$), for Péclet numbers ${\text {Pe}}=0.5$, 2.0 and 10.0

The Péclet number ${\text {Pe}} = (v{\varDelta }x) / D$ classifies the system as diffusion-dominated (${\text {Pe}} < 1$) or advection-dominated (${\text {Pe}} \gg 1$). For $D = 10^{-5}~{\hbox {m}}^2/\hbox {s}$ and ${\varDelta }x =0.01~\hbox {m}$, we vary the velocity v to get ${\text {Pe}} = 0.5$, 2.0 and 10.0. Figure 8 provides a comparison between the deterministic and stochastic models. The latter involves $M =10{,}000$ particles; and the number of particles in each mesh element is normalized by M. Good agreement is observed for a range of Péclet number.

The discrepancy error $\epsilon$ provides a quantitative comparison between the two models. It is defined as the $\ell ^2$ norm of the difference between their respective solutions, $\epsilon = [ \sum _{i=1}^{N_x}(u_i-U_i)^2 ]^{1/2}$. Here $u_i$ is the deterministic solution at $x=(i-0.5){\varDelta }x$ and $U_i$ is the normalized number of particles within the ith mesh element, $(i-1){\varDelta }x \le x \le i{\varDelta }x$. The stochastic simulation is not affected by the mesh size ${\varDelta }x$; the mesh is only used to count the number of particles in each mesh element so that visual comparison can be made. Figure 8 exhibits $\epsilon$ as function of the number of particles M, where $M = 2^p\times 1000$ with $p =1,\ldots ,8$, for various Péclet numbers. As expected, the error decreases with M in all Pe regimes.

1.2 One-dimensional simulations

We consider a one-dimensional analog of the Krogh tissue cylinder model (Krogh 1922), in which the tissue is perfused by parallel blood vessels. According to this model, the tissue is composed of cylindrical regions (of radius R), such is that each region is fed by a single blood vessel and effectively isolated from the adjacent regions by a cylindrical surface of symmetry. Without loss of generality, we replace the cylindrical domain with a one-dimensional geometry shown in Fig. 9.

Fig. 9

One-dimensional simulation domain motivated by the Krogh tissue cylinder model

Fig. 10

A leukocyte chasing a bacterium in one spatial dimension. Left column: the deterministic (solid lines) and hybrid-method (diamonds) solutions at times $t = 50~\hbox {s}$ and 100 s. Right column: a the chemoattractant concentration at capture time, b representative trajectories of the leukocytes and bacteria up to the capture time, and c histogram of the capture times from 100 hybrid one-dimensional simulations

As a result, significant mathematical simplifications can be made in the analysis of the model (Keller and Segel 1971). The one-dimensional version of (5) is solved numerically on the spatial domain ${\varOmega }= [x_{\text {sk}}, x_{\text {cap}}]$ of length $R = x_{\text {cap}} - x_{\text {sk}}$, where $x_{\text {sk}} = 0$ and $x_{\text {cap}}$ are the locations of the skin and the capillary/venule, respectively.

In the one-dimensional version of the initial and boundary conditions from Sect. 3.5, we set $h_a =1$ and ${\mathcal {A}}_0 =10^{-5}~\hbox {M}$, with values of the remaining parameters listed in Table 1.

1.2.1 One leukocyte chasing one bacterium

In this numerical experiment, one bacterium is released at $x_{\text {sk}} = 0~\hbox {m}$ and one leukocyte at $x_{\text {cap}} = 10^{-4}~\hbox {m}$. The interval ${\varOmega }$ is divided into 200 mesh elements. This setting provides a simplified representation of a neutrophil chasing a bacterium in the lab (Rogers xxxx). The reaction-rate value in Table 1 is such that neither the leukocyte nor the bacterium grows and/or decays during the time interval $t \in [0,1000~{\text {s}}]$. We further assume that, upon encounter (numerically, once they enter the same mesh element), both the bacteria and leukocyte die immediately. To keep a sharper concentration gradient and save the simulation time, we set $D = 10^{-12}~{\hbox {m}}^2/\hbox {s}$ in this single test.

Fig. 11

One-dimensional deterministic (solid lines) and hybrid (diamonds) simulations of inflammation process with an intermediate-size population of leukocytes ($N_{{\text {L}}_0} = 8000$) and bacteria $N_{{\text {B}}_0} = 8000$. Left column: the solution at times $t = 1000~\hbox {s}$ and 2000 s. Right column: a change in the size of populations over simulation time, b zoomed changes during early times, and c trajectories of three randomly selected leukocytes

Fig. 12

One-dimensional deterministic (solid lines) and hybrid (diamonds) simulations of inflammation process with a large population of leukocytes ($N_{{\text {L}}_0} = 8\times 10^5$) and bacteria $N_{{\text {B}}_0} = 8\times 10^5$. Left column: the solution at times $t = 1000~\hbox {s}$ and 2000 s. Solid lines are the profile of deterministic solution and stars are the hybrid-method results. Right column: a change in the size of populations over simulation time, b zoomed-in changes during early times, and c trajectories of three randomly selected leukocytes

Left column of Fig. 10 reveals that the deterministic PDE solution does not provide an accurate description of this experiment. Our hybrid method, on the other hand, accurately captures the microscopic random dynamics of the cells. Right column of Fig. 10 shows one trajectory of the leukocyte “chasing” the bacterium. The Brownian motion of the leukocyte and the bacterium introduces strong stochastic effects in the trajectory, which is naturally absent in the deterministic solution. This is an extreme scenario where different scales invoke the need for different simulation approaches. Concentration, as an average (macroscopic) description generally used for large number of particles, becomes an inadequate quantity to represent the dynamic of single particles (bacterium and leukocyte).

The last frame of Fig. 10 shows the histogram of capture time from 100 hybrid simulations. The average capture time of these 100 simulations is 351.3035 s. Since the leukocyte has less space to move about in one dimension than in two, it captures the bacterium faster.

1.2.2 Large populations of leukocytes and bacteria

We examine the situation where relatively large numbers of bacteria and leukocytes are involved in the inflammation process. The latter takes place in the domain ${\varOmega }= [0,10^{-3}~{\text {m}}]$, which is divided into 200 mesh elements. Populations of intermediate size consist of $N_{{\text {L}}_0} = 8000$ leukocytes and $N_{{\text {B}}_0} = 8000$ bacteria. Figure 11 shows two profiles of the chemoattractant concentration and the normalized numbers of bacteria and leukocytes. Predictions of the two methods are similar, except for some stochastic fluctuations in the numbers of cells. Furthermore, the hybrid method yields more detailed information about the reactions between leukocytes and bacteria than its PDE-based counterpart. The reduction in the population size with time is reported in Fig. 11a, b. During the early stages of inflammation, before leukocytes reach bacteria near the wound, the population-size changes are mostly due to the growth of bacteria and the natural death of leukocytes. Three representative leukocyte trajectories are shown in Fig. 11c.

Fig. 13

Two-dimensional simulations of chemotaxis–motility–reactions of 5 particles. Top row: the chemoattractant concentration, ${\mathcal {A}}({\mathbf {x}},t)$, at $t=1000~\hbox {s}$. Bottom row: trajectories of the leukocytes and bacteria at times $t = 500~\hbox {s}$ and 1000 s, respectively. The lighter colors represent positions at earlier times, and the darker colors represent positions at later times. Red diamonds are leukocytes released from $x_{{\text {L}}_1}$, green diamonds are from $x_{{\text {L}}_2}$, blue diamonds are from $x_{{\text {L}}_3}$, and cyan diamonds represent bacteria traces

We repeat this experiment but for large populations of leukocytes and bacteria, $N_{{\text {L}}_0} = N_{{\text {B}}_0} = 8 \times 10^5$. As one would expect, the PDE-based and hybrid models now yield virtually indistinguishable results (Fig. 12). This proves our hybrid method to be accurate as well as flexible for multiscale simulations in multi-species biological systems. Moreover, reactions, population changes and individual cell behavior can be observed on the fly (Fig. 12).

Fig. 14

Two-dimensional simulations of chemotaxis–motility–reactions of $8\times 10^5$ particles. Temporal snapshots, at time $t=1000~\hbox {s}$, of a the chemoattractant concentration and the numbers of b leukocytes and c bacteria normalized with their initial population numbers

Fig. 15

Two-dimensional simulations of chemotaxis–motility–reactions of $8\times 10^5$ particles. Top row: trajectories of the leukocytes and bacteria at times $t = 500~\hbox {s}$ and 1000 s, respectively. The lighter colors represent positions at earlier times, and the darker colors represent positions at later times. Red diamonds are leukocytes released from $x_{{\text {L}}_1}$, green diamonds are from $x_{{\text {L}}_2}$, blue diamonds are from $x_{{\text {L}}_3}$, and cyan diamonds represent bacteria traces. Bottom row: change in the size of populations over simulation time and zoomed changes during early times

1.3 Two-dimensional simulations

1.3.1 Small populations of leukocytes and bacteria

Leukocytes enter the tissue from three locations along the blood vessel, as shown in Fig. 3, with the fixed number of leukocytes at each site $N_{{\text {L}}_1} = 2$, $N_{{\text {L}}_2}=2$ and $N_{{\text {L}}_3}=1$, respectively. The total of $N_{{\text {B}}} = 5$ bacteria invade from the wound. We set $x_{\text {out}} = 4\times 10^{-4}~\hbox {m}$ and $y_{\text {cap}} = 10^{-4}~\hbox {m}$.

The simulation results are presented in Fig. 13 up to time $t = 1000~\hbox {s}$. The leukocytes disperse over the domain as time goes by, showing strong stochastic dynamics due to motility, but tend to migrate to the wound due to the chemotaxis effects.

1.3.2 Large populations of leukocytes and bacteria

The total of $N_{{\text {L}}_0} = M_B = 8\times 10^5$ leukocytes enter the tissue from the blood vessel: $N_{{\text {L}}_1} = 3\times 10^5$ at $x_{{\text {L}}_1}$, $N_{{\text {L}}_2} = 3\times 10^5$ at $x_{{\text {L}}_2}$ and $N_{{\text {L}}_3} = 2\times 10^5$ at $x_{{\text {L}}_3}$. The total of $N_{{\text {B}}_0} = 8\times 10^5$ invades from the wound. This setting gives rise to the collective behavior of leukocytes (Fig. 14), in which the random fluctuations observed in the previous section are smoothed out. The normalized particle-number map is less chaotic and smoother than that in Fig. 6. More leukocytes congregate in the left part of the simulation domain than on the right because more of them are released from $x_{{\text {L}}_1}$ than from $x_{{\text {L}}_3}$. This is well understood in terms of diffusion of concentration. The reaction zone is more pronounced in the results for large populations.

The individual behavior of leukocytes has the same pattern as before (Fig. 15). Changes in the population size have the same trends as those observed in the intermediate-size populations. However, the reactions are more intense during the same zoomed-in time interval due to a larger number of leukocytes.