Team:ETH Zurich/Modeling/Reaction-diffusion

"What I cannot create I do not understand."
- Richard Feynmann

Reaction-diffusion Models

Introduction

While single-cell models are useful for correctly implementing and debugging chemical reaction models, they are not sufficient to fully understand the real-life functionality of our system. Since an essential part of our system is increasing the perceived concentrations of lactate and AHL through co-localization, it is necessary to model the concentrations the chemical species though a reaction-diffusion system.

"Doughnut" model of our system: derivation

Figure 1: "Doughnut" model for our system.

Geometry: Abstraction of a population of discrete E. coli to a connected uniform space

The well in which our reaction takes place is a box of length 200μm. In the middle of the well, we have a circle representing our target cell. Surrounding this is a ring of width equal to the diameter of an E. coli cell representing the union of all E. coli cells bound to the target cell, highlighted in Figure 1, which we call the doughnut. We assume that the occupancy of the doughnut is constant and maximal (25 cells in a 2D model and 125 in a 3D model). Finally, the remaining space represents both the medium and free-floating E. coli cells, which we refer to as the bulk. All species are able to diffuse freely through this medium and rates for our reactions are adjusted based on the size of the E. coli population at a given time.

Logistic growth of E. coli

If we let \(n_\text{bulk}:[0,t_\text{sim}]\longrightarrow \mathbb N\) be the function giving the number of E. coli in the bulk at time \(t\), with an initial population of size \(n_0\), a carrying capacity of \(K\), and a growth rate of \(R\), then $$n_\text{bulk}(t) = \frac{n_0 K e^{Rt}}{K + n_0(e^{Rt}-1)}$$ so with an E. coli doubling time of \(t_2\) (which we assume is 30 minutes), we can solve for \(R\) $$R = \frac{\log2 + \log(K+n_0) - \log(K-2n_0)}{t_2}$$

Concentration correction for differing volumes

Since the areas of the doughnut and bulk are not equal to the true areas of the spaces they are representing, the concentrations of the chemical species in those areas need to be multiplied by a correcting factor to obtain their correct values. Let \(n\) represent the number of chemical species and \(t_\text{sim}\) be our total simulation time. Let \(\mathbf X:[0,t_\text{sim}]\longrightarrow \mathbb R^n\) be a function representing the molar concentrations of our chemical species over the simulation period. Then our system of non-linear ordinary differential equations (ODEs) can be represented by the following equation $$\left.\frac{d\mathbf X}{dt}\right|_{(\mathbf X(t),t)} = f(\mathbf X(t),t)$$ the units of which are \(\frac{\text{mol}}{L_\textit{E. coli}}\) for each component of \(\mathbf X\). For a more detailed description of the system, its assumptions and its derivation, please see the pages on the lactate module, AHL module and the combined model.

We can define a new function \(\mathbf Y(t) := \nu(t)\mathbf X(t)\) representing the concentrations of the species within our simulated bulk, where $$\nu_\text{bulk}(t) := \frac{n(t)A_\textit{E. coli}}{A_\text{bulk}}$$ $$\nu_\text{doughnut}(t) := \frac{n_\text{doughnut}A_\textit{E. coli}}{A_\text{doughnut}}$$ are concentration correction factors corresponding to the abstract spaces. Our original ODE system after this change of variables is then $$\left.\frac{d\frac{\mathbf Y}{\nu}}{dt}\right|_{(\mathbf Y(t),t)} = f\left(\frac{\mathbf Y(t)}{\nu(t)},t\right)$$ If we solve for \(\frac{d\mathbf Y}{dt}\), we get $$\left.\frac{d\mathbf Y}{dt}\right|_{(\mathbf Y(t),t)} = \nu(t)\left(f\left(\frac{\mathbf Y(t)}{\nu(t)},t\right) - \mathbf Y\frac{d}{dt}(\nu(t))^{-1}\right)$$ We can then derive ODE models for the doughnut and the bulk by plugging in their respective concentration correction factors. \begin{align*} \left.\frac{d\mathbf Y_\text{bulk}}{dt}\right|_{(\mathbf Y(t),t)} &= \nu_\text{bulk}(t)\left(f\left(\frac{\mathbf Y(t)}{\nu_\text{bulk}(t)},t\right) - R\left(\frac{n(t)}{K}-1\right)\mathbf Y(t)\right)\\ \left.\frac{d\mathbf Y_\text{doughnut}}{dt}\right|_{(\mathbf Y(t),t)} &= \frac{n_\text{bulk}A_\textit{E. coli}}{A_\text{doughnut}}\left(f\left(\frac{A_\text{doughnut}}{n_\text{doughnut}A_\textit{E. coli}}\mathbf Y(t),t\right)\right) \end{align*}

Diffusion and transport of chemical species

Under alkaline conditions, E. coli actively import lactate via a proton-motive symporter lldP.[] Thus, a cross-membrane transport reaction had to be implemented. Since this is not possible directly in COMSOL, we had to model lactate in two states. Suppose our reference is the subspace of the interiors of the E. coli. We then defined the two states \(Lact_\text{in}\) and \(Lact_\text{out}\), denoting intracellular and extracellular lactate, respectively. \(Lact_\text{out}\) is produced by the target cell and can diffuse freely though the medium and all membranes. \(Lact_\text{in}\) is in equilibrium with \(Lact_\text{out}\) with rate constants set to maintain a 20-fold difference of lactate concentration between interior and exterior. $$ Lact_\text{out} \mathop{\mathop{\xrightarrow{\hspace{4em}}}^{\xleftarrow{\hspace{4em}}}}_{k_\text{in}}^{k_{\mathrm{ext}}} Lact_\text{in} \qquad \frac{k_\text{in}}{k_\text{out}}\approx 20 $$ In addition, only the \(Lact_\text{in}\) state can react with the other chemical species in the E. coli.

AHL is able to freely diffusion in the medium and across membranes. The effective diffusion coefficients of AHL and Lactate through the E. coli membrane \(D_e\) were approximated by the method proposed by [Stewart 2003] as a fraction of their respective diffusion coefficients in water \(D_{aq}\) by the relation $$\frac{D_e}{D_{aq}}\approx 0.25$$ For all other chemical species, diffusion across the doughnut membrane represents the unbinding of a cell from the target and its introduction to the bulk. Due to the very high affinity of annexin V to phosphatidylserine (\(K_m = \)[]), we can make the assumption that another cell will bind in its place immediately, and thus, keep the doughnut population constant. The diffusion coefficients across the doughnut membrane are determined by the \(k_\text{on}\) and \(k_\text{off}\) rates of the annexin V-phosphatidylserine complex. Since diffusion is concentration-dependent, and the concentrations in the bulk need to be corrected, a correction must also be applied to the diffusion coefficients. Since direction-dependent diffusion coefficients are not supported by our multiphysics software, an approach similar to that for lactate was used, with the only difference being that since these species are only located in cell cytoplasm, the concentration correction function \(\nu\) has to be applied. Thus, for each component \(X^i\), the following equilibrium will be simulated $$ X^i_\text{out} \mathop{\mathop{\xrightarrow{\hspace{4em}}}^{\xleftarrow{\hspace{4em}}}}_{k_\text{on}}^{k_{\mathrm{off}}} X^i_\text{in} $$

Testing our system's logic through controls

Figure 2: Test well with four control wells.

To test whether our system acts as an AND gate on our two inputs (higher lactate production and co-localization signals), we combinatorially tested our system in environments with high vs. low lactate production and E. coli co-localization vs. dispersion, along with a fifth well with no target cell present.

Results

Parameter search on a simplified reaction-diffusion model