Learning black- and gray-box chemotactic PDEs/closures from agent based Monte Carlo simulation data

Abstract

We propose a machine learning framework for the data-driven discovery of macroscopic chemotactic Partial Differential Equations (PDEs)—and the closures that lead to them- from high-fidelity, individual-based stochastic simulations of Escherichia coli bacterial motility. The fine scale, chemomechanical, hybrid (continuum—Monte Carlo) simulation model embodies the underlying biophysics, and its parameters are informed from experimental observations of individual cells. Using a parsimonious set of collective observables, we learn effective, coarse-grained “Keller–Segel class” chemotactic PDEs using machine learning regressors: (a) (shallow) feedforward neural networks and (b) Gaussian Processes. The learned laws can be black-box (when no prior knowledge about the PDE law structure is assumed) or gray-box when parts of the equation (e.g. the pure diffusion part) is known and “hardwired” in the regression process. More importantly, we discuss data-driven corrections (both additive and functional), to analytically known, approximate closures.

Appendices

Appendix A: Details on the Monte Carlo chemotaxis model

Our microscopic, agent-based model is based on the work of Othmer and Schaap (1998); Othmer et al. (2013). Each bacterium is modelled as having six flagellae; special care has been taken in modelling the direction of the rotation of the flagellar filaments, as this constitutes the basis of chemotaxis (Larsen et al. 1974; Spiro et al. 1997). Following Scharf et al. (1998), the motor dynamics are described by a two-state system modelling the transition rates (transition probabilities per unit time) between CCW and CW (counter-clockwise and clockwise, respectively) rotation for each flagellum. These are characterized by an exponential distribution of time intervals in each state (Turner et al. 1996). Let us denote by $k^+$ ($k^-$) the transition rate from CCW to CW (CW to CCW). Then, the bias of CW, i.e. the fraction of time that a flagellum rotates CW is $p_{CW}=\frac{k^+}{k^{+}+k^-}$, $p_{CCW}=1-p_{CW}$. The reversal frequency in the direction of rotation of the flagellar motors is Turner et al. (1996):

$$\begin{aligned} \rho =p_{CCW} k_{+} + p_{CW} k_{-} = \frac{2k^+k^-}{k^++k^-}. \end{aligned}$$

(A.1)

and the rate constants are given by:

$$\begin{aligned} k_{+}=\rho /(2 p_{CCW}), k_{-}=\rho /(2 p_{CW}) \end{aligned}$$

(A.2)

For a cell with N flagellae, the total CCW bias of the cell is given by Spiro et al. (1997):

$$\begin{aligned} P_{CCW} = \sum _{j=\theta }^{N} \left( {\begin{array}{c}N\\ j\end{array}}\right) p_{CCW}^{j}(1-p_{CCW})^{N-j}. \end{aligned}$$

(A.3)

For $N=6$, $p_{CCW}=0.64$ and $\theta =N/2$, we get $P_{CCW}\sim 0.87$ suggesting that the cell spends around 90% of the time running (Spiro et al. 1997). This result is in line with experimental observations for the motility of wild-type E. coli in the absence of changes of the substrate, where the mean run (swimming) periods are $\sim 1$ s and the tumble periods $\sim 0.1$ s (for the strain AW405 in dilute phosphate buffer at $32^{\circ }\hbox {C}$) (Ishihara et al. 1983).

Experimental studies have shown that these rates depend on the CheY-P concentration, say C. In particular, Cluzel et al. (2000) have shown that the dependence of CW bias (between the values 0.1 and 0.9) to C can be approximated by a Hill function with a coefficient $H \sim 10.3 \pm 1.1$, with a dissociation constant $K_d= 3.1$ mM/s. Thus, the CW bias reads:

$$\begin{aligned} p_{CW} = \frac{C^H}{K_{d}^H +C^H}. \end{aligned}$$

(A.4)

Based on the above findings, the transition rates $k^+$, $k^-$ are given by Setayeshgar et al. (2005):

$$\begin{aligned} k^+= & {} \frac{H C^{H-1}}{K_{d}^H + C^H},\end{aligned}$$

(A.5)

$$\begin{aligned} k^-= & {} \frac{1}{C} \frac{H K_{d}^{H}}{K_{d}^H + C^H}. \end{aligned}$$

(A.6)

Thus, based on the model formulation and nominal values of the parameters, the expected fraction of time spent in the CCW state in the absence of stimulus for each cell from kMC simulations is $\sim$ 0.855, close enough to the one observed experimentally. In the absence of spatial variations in the chemoattractant (or repellent) profile, the rotation of the flagellar filament is biased towards the CCW direction (that is, the probability of CCW rotation of a flagellum is higher than that of CW rotation), when viewed along the helix axis towards the point of insertion in the cell (Larsen et al. 1974). This bias depends on the type of bacterial strain and the temperature; for the wild-type strain AW405, it has been found that the average value of the CCW bias is 0.64 at $32^{\circ }\hbox {C}$ (Larsen et al. 1974; Block et al. 1982). When the majority of the flagellar filaments rotate CCW (CW) the cell swims (tumbles).

Monte Carlo Model

Appendix B: Determination of the parameters of the macroscopic PDE for bacterial density evolution

1.1 Appendix B.1: Determination of the diffusion coefficient

An estimation of the diffusion coefficient for cell motility in the absence of stimulus can be attempted following two paths. From a microscopic point of view, considering a random walk simulation, the mean free path i.e., the swimming distance without any change in the direction is given by $\delta r = \tau \cdot \bar{v}$, where $\tau$ is the mean time of swimming in one direction. Considering n such time steps in time t (i.e. $n=t/\tau$), the total mean-squared displacement $\Delta r (t)^2$ at a certain time (t) is given by the Einstein relation (Berg and Turner 1990):

$$\begin{aligned} \langle \Delta r^2(t)\rangle =2 D_m t \approx 2 n \delta r^2 = 2 \tau \bar{v}^2 t, \end{aligned}$$

(B.1)

which is valid for $t>> \tau$, where, $\tau$ is the characteristic time scale. Here, the value of $D_m$ is estimated from our Monte Carlo simulations in the absence of stimulus (we have set $s(x)=1$, $\forall x$), by tracking the trajectories of 1000 cells for a time period of $2000~\textrm{s}$. The cells are initially positioned at the middle of the domain, all initialized at the tumbling phase, with $u_1(0)=0$ (no excitation), and adapted with $u_2(0)=f(s(x))$, with a constant velocity of $\bar{v}=0.003~\hbox {cm}/\hbox {s}$ (as in Berg and Turner 1990). Figure 5 depicts the average of the square distance as a function of time. By least-squares, we get $\hat{D}_m \approx 9 \cdot 10^{-6}~\hbox {cm}^{2}/\hbox {s}$. This is in good agreement with experimental observations for the E. coli motility (see Berg and Brown 1972; Berg and Turner 1990; Spiro et al. 1997; Cluzel et al. 2000).

Fig. 5

The average square displacement (over all 1000 cells) as a function of time in the absence of stimulus (gradient of chemoattractant), when all cells are initialized in the tumbling phase with $u_1(0)=0$, i.e. without excitation and fully adapted, i.e. $u_2(0)=f(s)$, and $\bar{v}=0.003~\hbox {cm}/\hbox {s}$

From a macroscopic point of view, one can estimate the diffusion coefficient $D_M$ from a linear curve fitting between $\frac{\partial b}{\partial t}$ and $\frac{\partial ^2 b}{\partial x^2}$ with finite difference approximations of temporal and spatial derivatives at the coarse-scale. Thus, by fixing a spatial gradient of chemo-nutrient profile to zero ($\nabla c = 0$), we can consider a simple diffusion equation with a constant diffusion coefficient, D:

$$\begin{aligned} \frac{\partial b}{\partial t} = D \nabla ^2 b. \end{aligned}$$

(B.2)

Finally, we note that the Einstein relation for the diffusion coefficient given by Eq. (B.1) can be approximated on average over a run and tumble period as:

$$\begin{aligned} <\Delta r^2>\approx \bar{v}^2 \bar{T}_{run}^2 =2 \bar{D}_m (\bar{T}_{run}+\bar{T}_{tumb}), \end{aligned}$$

(B.3)

where ${\bar{T}_{run}, \bar{T}_{tumb}}$ denote the average duration of swimming and tumbling periods, respectively, and $\bar{v}$ is the average swimming speed. Thus, based on Eq. (B.3) and assuming that the tumbling duration is negligible compared to the swimming duration (as assumed for the derivation of the generalized Keller–Segel theory embodied in Eq. (3)), an approximation of the diffusion coefficient is given by:

$$\begin{aligned} \bar{D}_m=\frac{\bar{v}^2}{2\lambda _0}, \quad \lambda _0=\bar{T}_{run}^{-1}. \end{aligned}$$

(B.4)

Hence, setting $\bar{D}_m=\hat{D}_m\approx 9 \cdot 10^{-6}~\hbox {cm}^2/\hbox {s}$, $\lambda _0=1~\hbox {s}^{-1}$ (in agreement with experimental findings), we get $\bar{v}=\sqrt{2}v=3 \sqrt{2} \cdot 9 \cdot 10^{-3}~\hbox {cm}/\hbox {s}$ as the average velocity appearing in Eq. (3).

1.2 Appendix B.2: Determination of the parameter c of the macroscopic PDE

As stated in Sect. 2, one of the assumptions for the derivation of the closed-form Keller–Segel Eq. (3) is the linear relation between the turning frequency $\lambda$, and the basal frequency $\lambda _0 \sim 1~\hbox {s}^{-1}$, i.e. for each cell at position x at time t, we have (see Eq. (5)):

$$\begin{aligned} \lambda (x,t)=\lambda _0-c u_1(x,t). \end{aligned}$$

(B.5)

For initial values $u_1(x,0)$, $u_2(x,0)$ for all cells (i.e. $\forall x \in \mathrm I\!R$), the analytical solution of the cartoon model (Eq. 2) is given by:

$$\begin{aligned} u_1(x,t)= & {} \frac{{\textrm{e}}^{-\frac{t}{\tau _{e}}}\,\left( {K_{s}}^2\,\tau _{a}-{K_{s}}^2\,\tau _{e} +s^2\,\tau _{a}-s^2\,\tau _{e}+2\,K_{s}\,s\,\tau _{a}-2\,K_{s}\,s\,\tau _{e}\right) }{{\left( K_{s}+s\right) }^2\,\left( \tau _{a}-\tau _{e}\right) }\,u_{1}(0,x) \nonumber \\{} & {} +\frac{{\textrm{e}}^{-\frac{t}{\tau _{e}}}\,(s^2\,\tau _{a}\,u_{2}-k\,s\,\tau _{a}+{K_{s}}^2\, \tau _{a}\,u_{2}(0,x)+2\,K_{s}\,s\,\tau _{a}\,u_{2}(0,x))+k\,s\,\tau _{a}\,}{{\left( K_{s}+s\right) }^2\, \left( \tau _{a}-\tau _{e}\right) }\nonumber \\{} & {} -\frac{\tau _{a}\,{\textrm{e}}^{-\frac{t}{\tau _{a}}}\,({K_{s}}^2\,u_{2}(0,x)-k\,s +s^2\,u_{2}(0,x)+2\,K_{s}\,s\,u_{2})+ k\,s\,\tau _a}{{\left( K_{s}+s\right) }^2\,\left( \tau _{a}-\tau _{e}\right) }, \end{aligned}$$

(B.6)

$$\begin{aligned} u_2(x,t)= & {} \frac{{\textrm{e}}^{-\frac{t}{\tau _{a}}}\,}{{\left( K_{s}+s\right) }}\,u_{2}(0,x) -\frac{k\,s({\textrm{e}}^{-\frac{t}{\tau _{a}}}-1)}{{\left( K_{s}+s\right) }^2}. \end{aligned}$$

(B.7)

Note that, if one sets as initial value $u_2(x,0)=f(s)$, then the second equation of the cartoon model (see Eq. (2)) gives $u_2(x,t)=f(s)$, $\forall t$ and the analytical solution for $u_1(x,t)$ is reduced to:

$$\begin{aligned} u_1(x,t)=u_{1}(x,0){\textrm{e}}^{-t/\tau _a}. \end{aligned}$$

(B.8)

To this end, the parameter c in Eq. (5) appearing in the Keller–Segel-class PDE given by Eq. (3) can be found with the aid of Monte Carlo simulations, by fixing $u_{1}(x,t)$, $\forall t$ to different relatively small values, say $u_1$ $\forall x$, and measuring, the number of turning events $\lambda (u_1)$; then the value of the parameter c can be estimated by least-squares. A different way would be to set an initial value for $u_1(x,0)$ (setting also as initial value $u_2(x,0)=f(s)$), run the Monte Carlo simulator, measure the turning frequencies $\lambda (u_1(t))$ and based on the above, the value of c can be again estimated with least-squares.

Here, to estimate c, we have fixed $u_1$ to the following values: $-0.02$, $-0.015$, $-0.01$, $-0.005$, 0, 0.005, 0.01, 0.015, 0.02, where the linear relation between $\lambda$ and $\lambda _0$ is valid, and we computed $\lambda (u_1)$ based on Monte Carlo simulations with 1000 cells for a time period of 2000s. For these values, $\hat{\lambda _0}\sim 1 s^{-1}$ (in a good agreement with the experimental findings) and $\hat{c}\sim 19.5$ (99% CI 18–21). We note that this value is consistent with what has been reported in other studies (Xue 2015). For our simulations with the macroscopic PDE, we have set $\hat{c}=20$.

Appendix C: Numerical details

See the Tables 2, 3 and 4.

Table 2 Relative % errors, for the reproduction case (corresponding to Fig. 3)

Table 3 Relative % error, for the testing case (corresponding to Fig. 4)

Table 4 ARD weights