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.

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.

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.

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}

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.

Model assumptions and limitations

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

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.