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).