Model

Introduction

In order to degrade drugs in hospital wastewater in a sustainable way, our biological system split into two distinct populations with different biological properties. Specifically the growth of the "stem" cell population may be faster that the degrader population. Here, degrader cells produce two enzymes, folC and CPG2, which can be lethal in high concentration.

We analyse the population dynamics for several scenarios. It could happen in hospital effluents, where drugs release depend of many parameters like human errors, schedule of health services, etc. Finally, we prove the effectiveness of our system by proving it can degrade drugs in many scenarios.

Model construction.

A first approach : the exponential phase

Let's denote $n_1(t)$ the population of degradation cells, $n_2(t)$ the population of "stem" cells and $m(t)$ the quantity of methotrexate in solution.

According to our experiments, and general theoretical consideration, we can write the following reactions : \begin{{ '{' }}array}{{ '{' }}clll} n_1 & \xrightarrow{{ '{' }}r_1} & n_1 & \text{{ '{' }}(degrader cells growth)} \\ n_2 & \xrightarrow{{ '{' }}r_2} & n_2 & \text{{ '{' }}("stem" cells growth.)} \\ n_1 + m & \xrightarrow{{ '{' }}-A} & n_1 & \text{{ '{' }}(degradation of methotrexate)} \\ \star & \xrightarrow{{ '{' }}p} & m & \text{{ '{' }}(methotrexate input)} \\ \end{{ '{' }}array}
The letter above the arrow is the rate at which the transformation happens. Here, all theses rates are assumed to be positive.

Now, in order to transform those reactions into mathematical solutions we use the law of mass action . This law states that transformations are proportional to the input reactants.
Therefore : \begin{{ '{' }}eqnarray} \frac{{ '{' }}dn_1}{{ '{' }}dt} & = & r_1 n_1 & (1.1)\\ \frac{{ '{' }}dn_2}{{ '{' }}dt} & = & r_2 n_2 & (1.2)\\ \frac{{ '{' }}dm}{{ '{' }}dt} & = & - A m n_1 + p & (1.3). \end{{ '{' }}eqnarray}

The environment carrying capacity

Over time and without methotrexate, both cell populations reach a stationary phase due to the finite carrying capacity $K$ of the environment. This new interaction explains that infinite population growth is unsustainable because the amount of renewable ressources in the environment is finite.

To this end, we multiply growth terms in equation (1.1) and (1.2) by $[1 - (n_1 + n_2)/K]$. This operation is known as "The Verhulst equation" .

However at the steady state, this new term prevent cells from switching from $n_1$ to $n_2$ or from $n_2$ to $n_1$. Therefore, we add switching terms in our equations.

So, we have new biological reactions : \begin{{ '{' }}array}{{ '{' }}clll} \star & \xrightarrow{{ '{' }}\text{{ '{' }}environment}} & n_1 & \text{{ '{' }}(carrying capacity of the environment)} \\ \star & \xrightarrow{{ '{' }}\text{{ '{' }}environment}} & n_2 & \text{{ '{' }}(carrying capacity of the environment)} \\ n_1 & \xrightarrow{{ '{' }}a} & n_2 & \text{{ '{' }}(switching from $n_1$ to $n_2$ at the steady state)} \\ n_2 & \xrightarrow{{ '{' }}b} & n_1 & \text{{ '{' }}(switching from $n_2$ to $n_1$ at the steady state)}. \\ \end{{ '{' }}array}
The letter above the arrow is the rate at which the transformation happens. Here, all theses rates are assumed to be positive.

Through the law of mass action, we transform these reactions into mathematical equations : \begin{{ '{' }}eqnarray} \frac{{ '{' }}dn_1}{{ '{' }}dt} & = & r_1 n_1 ( 1 - \frac{{ '{' }}n_1 + n_2}{{ '{' }}K}) + b n_2 - a n_1 & (eqn.A)\\ \frac{{ '{' }}dn_2}{{ '{' }}dt} & = & r_2 n_2 ( 1 - \frac{{ '{' }}n_1 + n_2}{{ '{' }}K}) + a n_1 - b n_2 & (eqn.B)\\ \frac{{ '{' }}dm}{{ '{' }}dt} & = & - A m n_1 + p. & (eqn.C) \end{{ '{' }}eqnarray}

Model assumptions and limitations

Our model relies on a few approximations and assumptions, which we found both necessary to keep the problem tractable, and realistic enough to attain our goal of giving us actionnable insight. We assume that all relevant reactions can be accounted for accurately with a small number of linear ordinary differential equations (ODE). The full pathway is of course more complex but we took the decision to abstract and simplify anytime we could afford to without compromising predictive power.

Furthermore, the linearity hypothesis only holds locally in the state space : as long as all parameters are "reasonnable" then the approximation error is largely negligible. However considering more extreme regimes would be both impractical (hard to gather reliable experimental data) and irrelevant in pratice for our applied and industrial setting.

While we have worked throughout our projet on a variety of models, we have found the ones presented here to be the most balanced in regards to the complexity / accuracy tradeoff, and most useful for effective decision making. In particular we expanded great effort in keeping our models analytically tractable where possible, as it is our opinion that analytical (ie. explicit) solutions are often much more enlightening to practitionners for further biological system design than numerical solutions.

Mathematical Analysis

The steady state

Figure 1 : Steady State without Methotrexate, a = 0.7, b = 0.3

Asymptotically, the population reaches steady state value $n_1 + n_2 = K$. Then the ratio between the "stem" cells and the degradation cells only depends on the switching parameters $a$ and $b$ (strain and pathway specific). $$ f = \frac{{ '{' }}n_1}{{ '{' }}n_2} = \frac{{ '{' }}b}{{ '{' }}a}. $$ Figure 1 illustrates how $n_1$ and $n_2$ eventually reach a steady state when long times are considered.
In the end, $$ n_1 = \frac{{ '{' }}fK}{{ '{' }}f + 1} $$ Moreover, this steady state is stable. It means that small pertubations around the equilibrium will not change the state of our system.

When the system reaches steady state, (eqn.C) becomes : $$ 0 = \frac{{ '{' }}dm}{{ '{' }}dt} = - A m n_1 + p = -A m \frac{{ '{' }}fK}{{ '{' }}f + 1} + p$$ and $$ m = \frac{{ '{' }}f + 1}{{ '{' }}fK} p A $$ Then, Methotrexate concentration converges to $$\frac{{ '{' }}f + 1}{{ '{' }}fK} p A$$

Proof

1) At the steady state $n_1 + n_2$ = K. So we can simplify (eqn.A) and (eqn.B) : \begin{{ '{' }}eqnarray} 0 & = & b n_2 - a n_1 \\ 0 & = & a n_1 - b n_2\\ \end{{ '{' }}eqnarray}

Therefore $$ b n_2 = a n_1 $$ and $$ \frac{{ '{' }}n_1}{{ '{' }}n_2} = \frac{{ '{' }}b}{{ '{' }}a} $$
2) At the steady state we have $$n_1 + n_2 = K \text{{ '{' }} and } f = \frac{{ '{' }}n_1}{{ '{' }}n_2} = \frac{{ '{' }}b}{{ '{' }}a}$$ Therefore $$\frac{{ '{' }}n_1}{{ '{' }}n_1}+ \frac{{ '{' }}n_2}{{ '{' }}n_1} = \frac{{ '{' }}K}{{ '{' }}n_1} $$ So $$1 + \frac{{ '{' }}1}{{ '{' }}f} = \frac{{ '{' }}K}{{ '{' }}n_1} $$ So $$ n_1 = \frac{{ '{' }}fK}{{ '{' }}f + 1}. $$ 3) Moreover this steady state is stable. In order to show that, let's take $\varepsilon > 0$ where $\varepsilon$ is a little perturbation. Then : \begin{{ '{' }}eqnarray} \frac{{ '{' }}dn_1}{{ '{' }}dt} & = & r_1n_1 (1 - \frac{{ '{' }}n_1 + n_2 + \varepsilon}{{ '{' }}K}) - an_1 + bn_2 \\ & = & -r_1n_1 \varepsilon + r_1n_1 (1 - \frac{{ '{' }}n_1 + n_2}{{ '{' }}K}) - an_1 + bn_2. \end{{ '{' }}eqnarray} Near the equilibrium : $$ n_1 + n_2 \approx K \text{{ '{' }} and } f = \frac{{ '{' }}n_1}{{ '{' }}n_2} \approx \frac{{ '{' }}b}{{ '{' }}a}.$$ Therefore : $$ \frac{{ '{' }}dn_1}{{ '{' }}dt} \approx -r_1n_1 \varepsilon < 0 $$ We see that the value is negative, so we conclude that the steady state is stable.

Towards a system to degrade more drugs

From a proof of concept to the general problem

Until now, our system was designed to only degradate methotrexate. But for us, methotrexate is just a proof of concept, a way to show we can degrade any drug with our system.
This is where our model becomes important because it can be adapted for any situation. This is the main motivation : be able to iterate with a predictive model driven approach faster than would be possible with only wetlab experiments.

Here again let's denote $n_1(t)$ the population of degradation cells, $n_2(t)$ the population of "stem" cells and $d(t)$ the quantity of targeted drug.
The main difference with our precedent model is that hazardous drugs may kill some of our cells. To represente that we indroduce two growth term. The first one $\mu^g$ when there are no drugs, and the second one $\mu^s$ when there are drugs in the environnement.
Specifically $\mu^s < \mu^g$.

The transformation can be represented by biological reactions : \begin{{ '{' }}array}{{ '{' }}clll} n_1 & \xrightarrow{{ '{' }}\mu^g_1} & n_1 & \text{{ '{' }}(degradation cells growth without drugs)} \\ n_2 & \xrightarrow{{ '{' }}\mu^g_2} & n_2 & \text{{ '{' }}("stem" cells growth without drugs.)} \\ n_1 & \xrightarrow{{ '{' }}\mu^s_1} & n_1 & \text{{ '{' }}(degradation cells growth with drugs)} \\ n_2 & \xrightarrow{{ '{' }}\mu^s_2} & n_2 & \text{{ '{' }}("stem" cells growth with drugs.)} \\ n_1 + m & \xrightarrow{{ '{' }}-A} & n_1 & \text{{ '{' }}(degradation of drugs)} \\ \star & \xrightarrow{{ '{' }}p} & m & \text{{ '{' }}(methotrexate input)} \\ \star & \xrightarrow{{ '{' }}\text{{ '{' }}environment}} & n_1 & \text{{ '{' }}(carrying capacity of the environment)} \\ \star & \xrightarrow{{ '{' }}\text{{ '{' }}environment}} & n_2 & \text{{ '{' }}(carrying capacity of the environment)} \\ n_1 & \xrightarrow{{ '{' }}a} & n_2 & \text{{ '{' }}(switching from $n_1$ to $n_2$ at the steady state)} \\ n_2 & \xrightarrow{{ '{' }}b} & n_1 & \text{{ '{' }}(switching from $n_2$ to $n_1$ at the steady state)}. \\ \end{{ '{' }}array}

We use the Law of Mass Action to translate theses reactions into mathematical equations : \begin{{ '{' }}eqnarray} \frac{{ '{' }}dn_1}{{ '{' }}dt} & = & \mu_1 n_1 ( 1 - \frac{{ '{' }}n_1 + n_2}{{ '{' }}K}) + b n_2 - a n_1 & (eqn.D)\\ \frac{{ '{' }}dn_2}{{ '{' }}dt} & = & \mu_2 n_2 ( 1 - \frac{{ '{' }}n_1 + n_2}{{ '{' }}K}) + a n_1 - b n_2 & (eqn.E)\\ \frac{{ '{' }}dm}{{ '{' }}dt} & = & - A m n_1 + p. & (eqn.F) \end{{ '{' }}eqnarray}

Response to an environmental shift

Figure 2 : The degradation of Methotrexate

Patra and Klumpp studied the dynamics of similar system when a population that has been growing without drugs is exposed to drugs. Before removing all kind of drugs (Fig.2).
There is a biphasis dynamic which translate into different phenotypes in the cell population. Authors said that :

"At time t~15 hours, the parameters were changed to those for stress condition. After the shift to stress conditions (by the addition of an antibiotic), the total population displays the biphasic decay behavior. In the fast-decaying phase, the decay of the total population is dominated by the death of normal cells, while in the second, slower-decaying phase, the total population consists predominantly of persister cells and the decay rate is governed by the death rate of the persisters.
The transition between the two different phases occurs when both subpopulation becomes equal in size" ( , p.5).

Figure from : Numerical integration of Eqs. (D,E) over a growth period of 15 hours (with μgn 2hr−1,μgn 0.2hr−1) followed by a stress period (with μgn −2hr−1,μgn −0.2hr−1) of another 15 hours. The switching rates were chosen to be a 0,001hr−1 and b 0,0001hr−1 . The killing curve of total population show two distinct phases, a fast-decaying phase and a slow-decaying phase. (b) Numerical integration of Eqs. (1) over a stress period of 15 hours followed by a regrowth period of another 15 hours. The regrowth curve of the total population shows two distinct phases, a slow-growing phase followed by a fast-growing phase. The parameters are the same as in (a). Here the "stem" cells are the normal cells, and the degradation cells are the persistence cells.

They even compute $T_s$, the time where transition occurs : $$T_s \approx \frac{{ '{' }}1}{{ '{' }}\Delta_s} \cdot ln(\frac{{ '{' }}n_2(0)}{{ '{' }}n_1(0)}) $$ where $\Delta_s = (\mu_1^s - a) - (\mu_2^s - b)$. We also note $\Delta_g = (\mu_2^g - b) - (\mu_1^g - a)$

Dynamics in periodically switching environment

In condition, drugs quantity depends of many factor : human errors, schedule of health services, etc. But even if that quantity isn't a constant, a looser but more realistic assumption would be that drug input is periodic.
Here, we will calculate the optimal parameters $a_{{ '{' }}opt}$ and $b_{{ '{' }}opt}$ in order to degrade methotrexate.

During one period, let's denote $t_s$ the duration of drugs presence in the environment and $t_g$ the duration without drugs. The average growth rates of our cells populations are : \begin{{ '{' }}eqnarray} <\mu_1> & = & \frac{{ '{' }}<\mu_1^g> t_g + <\mu_1^s> t_s}{{ '{' }}t_g + t_s} \\ <\mu_2> & = & \frac{{ '{' }}<\mu_2^g> t_g + <\mu_2^s> t_s}{{ '{' }}t_g + t_s} \\ \end{{ '{' }}eqnarray} where $<x>$ is the average growth rate of x.

Let's denote $\mu = \mu_1 + \mu_2$ the global growth rate. Then $$ <\mu> = \frac{{ '{' }}<\mu_1^g> t_g + <\mu_1^s> t_s}{{ '{' }}t_g + t_s} + \frac{{ '{' }}<\mu_2^g> t_g + <\mu_2^s> t_s}{{ '{' }}t_g + t_s} $$ Then, we use results in Patra and Klumpp [3] to compute that : $$ <\mu> = \frac{{ '{' }}(\mu_2^g - b) t_g + (\mu_1^s - a) t_s}{{ '{' }}t_g + t_s} + \frac{{ '{' }}ln(\frac{{ '{' }}ab}{{ '{' }}\Delta_s \Delta_g})}{{ '{' }}t_g + t_s} $$ We consider that $a$ and $b$ are optimal when they maximize $<\mu>$. That means that : $$ \frac{{ '{' }}1}{{ '{' }}a_{{ '{' }}opt}} = t_s - \frac{{ '{' }}1}{{ '{' }}\Delta_s} + \frac{{ '{' }}1}{{ '{' }}\Delta_g}$$ $$ \frac{{ '{' }}1}{{ '{' }}b_{{ '{' }}opt}} = t_g - \frac{{ '{' }}1}{{ '{' }}\Delta_g} + \frac{{ '{' }}1}{{ '{' }}\Delta_s} $$ Thus, we have found $a_{{ '{' }}opt}$ and $b_{{ '{' }}opt}$ in order to degrade drugs. If we had to build such a biological system, we would aim for these values when assessing possible promoters.

Team:GO Paris-Saclay/model