Reaction-diffusion Models

Overview

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 to our model.

The results from these simulations show that due to our experimental environment, the production of an AHL gradient oriented towards the outer edge of the well (i.e. higher concentration near the target cell) is not possible due to the fast diffusion of AHL. However, the greater capture of lactate by cells bound to the target relative to free-floating cells creates a concentration gradient of bound lactate and LuxR which translates to a difference in GFP signal between cells in the two locations. Combinatorially testing cases where E. coli cells bind or float around high and low lactate producers shows that our system acts as an AND gate of the two cancer markers, with E. coli bound to a high lactate producer fluorescing at orders of magnitude greater strength.

The "doughnut" model

Introduction

After implementing the single-cell models in the COMSOL multiphysics software package in an analogous geometry as a control, 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. To address this, the discrete circles were abstracted to connected spaces representing the union of a collection of cells 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 is referred to as the "doughnut" model. The AHL and lactate modules were implemented and tested separately before being combined into a full model.

All mathematical derivations are done in general, with values for variables being defined in the definition of our model geometry and assumptions, and in the parameters page.

Abstraction of discrete E. coli to a connected 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 surrounded by a ring (called the doughnut) of width equal to the diameter of an E. coli cell representing the union of all bound E. coli cells (see Figure 1). 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.

Assumptions

1. Uniform distribution of cells in their respective spaces
2. Cell density in the doughnut is constant and maximal
3. A cell unbinding from the target is immediately replaced
4. Logistic growth of E. coli in the bulk
5. E. coli doubling time of 30 minutes
6. One non-dividing target cell per well
7. Experiment ends when the size of the E. coli population reaches 95% of the carrying capacity
8. The carrying capacity of the well is the maximum number of cells that can physically fit in the well

AHL Module

Initial model: Constant LuxR, constant population size

The initial characterization of the system was done in a model in which the concentration of LuxR was constant and corresponded to the concentration generated by the lactate module in response to a high-lactate-producing target cell. In addition, the population sizes of the E. coli in the doughnut and bulk were kept constant at 125 and 1000 cells, respectively. Under this framework, the toggling of riboregulation of the LuxI promoter and AHL degradation by AiiA were two binary-valued parameters. The final ratio of GFP concentration between a cell in the doughnut and the bulk was used as a success measure.

• Without AiiA
• Error creating thumbnail: File missing

• With AiiA
• Error creating thumbnail: File missing

The results from this simulation indicate that the use of riboregulation is a sufficient and necessary condition for a difference in final GFP concentration to be produced. The presence of AHL degradation lowers the absolute values of the GFP concentrations, and lowers the ratio of final concentrations by almost an order of magnitude. It is also observed that, due to the small well size, the diffusion of AHL and loss of a concentration gradient occurs in less than 1 minute. These results indicate that the expression of AiiA does not generate a concentration gradient of AHL as originally planned.

In addition to testing the influence of riboregulation and AHL degradation, we also explored the effects of bulk population size on the final result. Starting from a size of 1000, a parameter sweep was performed in intervals of 500 cells until the E. coli carrying capacity of 12,000 cells was reached (see below). A linear decrease in the GFP concentration ratio was observed.

Lactate Module

Model Extensions

Logistic growth of E. coli

Let us define \(n_\text{bulk}:[0,t_\text{sim}]\longrightarrow \mathbb N\) as the function representing the size of the E. coli population in the bulk at time \(t\), with an initial population size of \(n_0\), a carrying capacity of \(K\), and a growth rate of \(R\). The closed-form equation of this function is then $$n_\text{bulk}(t) = \frac{n_0 K e^{Rt}}{K + n_0(e^{Rt}-1)}$$

If we denote the doubling time of the population by \(t_2\), we can solve for the growth rate \(R\) and define a new constant \(g\) s.t. $$R = \frac{\log2 + \log(K+n_0) - \log(K-2n_0)}{t_2} =: \frac{\log g}{t_2}$$ When we plug this value back into the original equation, it simplifies to $$n_\text{bulk}(t) = \frac{n_0 K g^{\frac{t}{t_2}}}{K + n_0(g^{\frac{t}{t_2}}-1)}$$ Finally, if we define the time we end our experiment, denoted \(t_\text{sim}\), as the time when the population reaches size \((1-\varepsilon)K\) for some small \(\varepsilon>0\), then we can solve for \(t_\text{sim}\) $$t_\text{sim} = \frac{1}{R}\log\frac{(K-n_0)(1-\varepsilon)}{n_0\varepsilon}$$

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 diffusion across barriers is implemented with equal diffusion coefficients for both directions in COMSOL, the transport of lactate across the doughnut membrane was implemented by modeling two different states of unbound lactate. If we let our reference space be the interior of the doughnut, we can define two pseudo-species \(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, while \(Lact_\text{in}\) is in equilibrium with \(Lact_\text{out}\) in the doughnut 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 $$ and cannot diffuse out to maintain its high concentration in the doughnut. Only the \(Lact_\text{in}\) state can react with the other chemical species in the doughnut and only \(Lact_\text{out}\) can react with other chemical species in the bulk.

Since AHL can freely diffuse across membranes, unbound AHL is only modeled as one state that can diffuse across all barriers. The effective diffusion coefficients of AHL and Lactate through the doughnut membrane \(D_e\) were approximated relative to their coefficients in water \(D_{aq}\) by the relation $$\frac{D_e}{D_{aq}}\approx 0.25$$ using the approximation proposed in a 2003 review by Stewart.

Correcting concentrations

Let \(\mathbf X:[0,t_\text{sim}]\longrightarrow \mathbb R^n\) be a function representing the molar concentrations of our chemical species in a real cell over the simulation period with units \(\frac{\text{mol}}{L}\) for each component of \(\mathbf X\). For a detailed description of the system, its assumptions and its derivation, please see the pages on the lactate module, AHL module and the combined model.

Since the areas of the doughnut and bulk are not equal to the areas of the spaces they are representing, a new ODE system equivalent to our model, but occurring in the model spaces, needs to be calculated. We define the function \(f\) such that our ODE system can be represented in the form \(\left.\frac{d\mathbf X}{dt}\right|_{(\mathbf X(t),t)} = f(\mathbf X(t),t)\). We can then define a new function \(\mathbf Y(t) := \nu(t)\mathbf X(t)\) representing the concentrations of the species within our model spaces, where $$\nu_\text{bulk}(t) := \frac{n(t)A_\textit{E. coli}}{A_\text{bulk}}\qquad \nu_\text{doughnut}(t) := \frac{n_\text{doughnut}A_\textit{E. coli}}{A_\text{doughnut}}$$ are concentration correction factors corresponding to our abstract spaces. Solving for \(\frac{d\mathbf Y}{dt}\) yields for each space yields \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*} Since \(AHL\) and \(Lact_\text{out}\) are the only species that are allowed to diffuse freely between the cells, we can obtain an equivalent system operating in the model spaces by multiplying the concentrations of AHL and \(Lact_\text{out}\) by \(\nu\) and dividing their derivatives by \(\nu\) to save on computing time.

Binding and unbinding of cells

None of the other chemical species are able to significantly diffuse across cell membranes in vivo. However, in order to model the unbinding and binding of cells to the target cell, an approach of implementing separate intracellular and extracellular pseudo-species was applied to the remaining chemical species. $$ X^i_\text{out} \mathop{\mathop{\xrightarrow{\hspace{4em}}}^{\xleftarrow{\hspace{4em}}}}_{k_\text{in}}^{k_{\mathrm{off}}} X^i_\text{in} $$ Unlike the model for lactate, the reaction rates for this equilibrium are governed not only by the concentrations of the relevant chemical species, but also by the number of cells that are part of the doughnut. So, a chemical analog of cell count has to be modelled. For this section, we are no longer assuming that the density of the doughnut is constant and maximal.

Binding of E. coli cells is through the binding of surface-exposed annexin V (A5) on their membrane to phosphatidylserine (PS) on the target cell surface to form the complex (APS). $$ A5 + PS \mathop{\mathop{\xrightarrow{\hspace{4em}}}^{\xleftarrow{\hspace{4em}}}}_{k_\text{on}}^{k_{\mathrm{off}}} APS $$ By our assumption that the target cell does not divide, the total amount of PS in the system is constant (\([PS]_t = [PS] + [APS]\)). If we let \([A5]_s\) denote the annexin V concentration per cell, the system can then be modeled by the ODE $$\frac{d[APS]}{dt} = k_\text{on}[A5]_s n_\text{bulk}(t)([PS]_t - [APS]) - k_\text{off}[APS]$$ If we assume that each E. coli cell has an equal amount of annexin V on its surface, then we can define the cell count as $$C(t) = \frac{[APS]}{[PS]_t}n_\text{around}$$ and an ODE system for the cell count $$\frac{dC}{dt} = k_\text{on}[A5]_s n_\text{bulk}(t)(n_\text{around} - C(t)) - k_\text{off}C(t)$$ Finally, we can use this to define an ODE system for the conversion of the intracellular and extracellular states of our chemical species $$\frac{dX^i_\text{out}}{dt} = k_\text{off}X^i_\text{in}C(t) - \alpha(t) k_\text{on}[A5]_sn_\text{bulk}(t)X^i_\text{out}(n_\text{around}-C(t))$$ where \(\alpha(t)\) is a rate correction term that has not yet been defined (will it ever be? who knows?).

Testing our system's logic

Figure 2: Test well with four control wells.

Ideally, our system should produce much stronger fluorescence only in the presence of both input signals (higher lactate production and co-localization signals). To test whether our system acts as an AND gate on these signals, we combinatorially checked all four cases of environments with high vs. low lactate production and E. coli co-localization vs. uniform distribution, 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, it was expected that a large time delay would exist between the high lactate and colocalization test and all others, or that a greater feedback would be produced in that environment.

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 small size of our well prevented this from occurring. Within a one minute period, any produced AHL quickly diffused 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 produced 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

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

Parameter search on a simplified reaction-diffusion model