"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 effective concentrations of lactate and AHL through co-localization, it is necessary to introduce a reaction-diffusion system.

After implementing the single-cell model in the COMSOL multiphysics software package in an analogous geometry, a more accurate geometry was constructed. Preliminary simulations were done in a geometry where a central circle representing the target cell was surrounded with smaller circles representing bound E. coli and circles uniformly distributed in the remaining space represented unbound E. coli. It was quickly realized, however, that this would not be an accurate representation of our test due to the rapid growth of the E. coli population. Thus, a geometry abstracting the discrete circles into continuous spaces was derived to allow for population growth to be simulated by modification of the ODE system over time. Owing to the shape of the layer of E. coli bound to the target cell, this geometry was referred to as the "doughnut" model.

The "doughnut" model

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

Figure 1: "Doughnut" model for our system.

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$$

Binding and unbinding of cells from the doughnut

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} $$ Each of the $X^i_\text{out}$ can diffuse freely through the doughnut membrane.

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. Although the positive feedback loop of the AHL module will eventually activate in all cases due to the leakiness of p_LuxR, the expectation is that this time of activation will be significantly reduced in environments with reduced lactate or no co-localization, thus acting as an AND gate of the two signals.

Results

Effects of colocalization

Figure 3: AHL concentration gradients for different production levels of the AHL degrader Aiia

Although our original plan was to strengthen quorum sensing between the doughnut cells through AHL degradation in the bulk, it became apparent that the high diffusion coefficient of AHL and the size of our well prevent this from occurring. Within a one minute period, any produced AHL quickly diffuses to produce a uniform distribution. The introduction of increased AHL degradation in the bulk led to an increased rate of AHL diffusion into the bulk and did not produce a gradient as desired.

Figure 4: LuxR concentration gradient

Fortunately, the co-localization of cells in the doughnut produce a substantial lactate and LuxR gradient due to lactate capture. Since p_LuxR is tied to both LuxI and GFP, this concentration gradient produces a stronger positive feedback response in these two proteins, thus producing the desired GFP concentration gradient.

Our system as an AND gate of two cancer markers

Figure 5:GFP concentrations in the test wells at time t_sim

Team:ETH Zurich/Modeling/Reaction-diffusion