Ratiometric control for differentiation of cell populations endowed with synthetic toggle switches

TL;DR

This paper develops and validates three control strategies to regulate the ratio of two cell populations with genetic toggle switches, considering practical constraints, using stochastic simulations.

Contribution

It introduces novel control methods for cell population ratios that incorporate realistic experimental constraints and validates them through in-silico simulations.

Findings

Successful regulation of cell population ratios in simulation.

Control strategies accommodate physical and technological constraints.

Validated effectiveness of controllers in stochastic agent-based models.

Abstract

We consider the problem of regulating by means of external control inputs the ratio of two cell populations. Specifically, we assume that these two cellular populations are composed of cells belonging to the same strain which embeds some bistable memory mechanism, e.g. a genetic toggle switch, allowing them to switch role from one population to another in response to some inputs. We present three control strategies to regulate the populations' ratio to arbitrary desired values which take also into account realistic physical and technological constraints occurring in experimental microfluidic platforms. The designed controllers are then validated in-silico using stochastic agent-based simulations.

Tables

Table 1. TABLE I: Value of the parameters of the cell population models (taken from [ 11 ] ).

${\kappa }_L^{\mathrm{m0}}$	$3.045\cdot {10}^{-1}$ mRNA ${\min }^{-1}$	${\gamma }_L^p$	$1.65\cdot {10}^{-2}$ ${\min }^{-1}$
${\kappa }_T^{\mathrm{m0}}$	$3.313\cdot {10}^{-1}$ mRNA ${\min }^{-1}$	${\gamma }_T^p$	$1.65\cdot {10}^{-2}$ ${\min }^{-1}$
${\kappa }_L^m$	$13.01$ mRNA ${\min }^{-1}$	${\theta }_{\mathrm{LacI}}$	$124.9$
${\kappa }_T^m$	$5.055$ mRNA ${\min }^{-1}$	${\eta }_{\mathrm{LacI}}$	$2.00$
${\kappa }_L^p$	$0.6606$ a.u. ${\text{mRNA}}^{-1}$ ${\min }^{-1}$	${\theta }_{\mathrm{TetR}}$	$76.40$
${\kappa }_T^p$	$0.5098$ a.u. ${\text{mRNA}}^{-1}$ ${\min }^{-1}$	${\eta }_{\mathrm{TetR}}$	$2.152$
$k_{\mathrm{aTc}}$	$4\cdot {10}^{-2}$ ${\min }^{-1}$	${\theta }_{\mathrm{aTc}}$	$35.98$
$k_{\mathrm{IPTG}}$	$4\cdot {10}^{-2}$ ${\min }^{-1}$	${\eta }_{\mathrm{aTc}}$	$2.00$
${\gamma }_L^m$	$1.386\cdot {10}^{-1}$ ${\min }^{-1}$	${\theta }_{\mathrm{IPTG}}$	$2.926\cdot {10}^{-1}$
${\gamma }_T^m$	$1.386\cdot {10}^{-1}$ ${\min }^{-1}$	${\eta }_{\mathrm{IPTG}}$	$2.00$

Table 2. TABLE II: Performance indices of the proposed controllers evaluated using Matlab (BSim, respectively) simulations.

Controller	$\overline{e}$	${\overline{e}}_f$	${\overline{t}}_s$  (min)
Bang-Bang	0.20 (0.13)	0.07 (0.06)	1077 (1185)
PI	0.17 (0.28)	0.02 (0.05)	563 (1020)
MPC	0.13 (0.22)	0.02 (0.05)	329 (820)

Equations

$$\frac{\mathrm{d}mRN{A}_{\mathrm{LacI}}^i}{\mathrm{dt}}={\kappa }_{\mathrm{L}}^{m0}+{\kappa }_{\mathrm{L}}^{\mathrm{m}}{\Phi }_{\mathrm{T}}(t)-{\gamma }_{\mathrm{L}}^{\mathrm{m}}mRN{A}_{\mathrm{LacI}}$$

$$\frac{\mathrm{d}mRN{A}_{\mathrm{TetR}}^i}{\mathrm{dt}}={\kappa }_{\mathrm{T}}^{m0}+{\kappa }_{\mathrm{T}}^{\mathrm{m}}{\Phi }_{\mathrm{L}}(t)-{\gamma }_{\mathrm{T}}^{\mathrm{m}}mRN{A}_{\mathrm{TetR}}$$

$$\frac{\mathrm{d}Lac{I}^i}{\mathrm{dt}}={\kappa }_{\mathrm{L}}^{\mathrm{p}}mRN{A}_{\mathrm{LacI}}^i-{\gamma }_{\mathrm{L}}^{\mathrm{p}}Lac{I}^i$$

$$\frac{\mathrm{d}Tet{R}^i}{\mathrm{dt}}={\kappa }_{\mathrm{T}}^{\mathrm{p}}mRN{A}_{\mathrm{TetR}}^i-{\gamma }_{\mathrm{T}}^{\mathrm{p}}Tet{R}^i$$

$$\frac{\mathrm{d}aT{c}^i}{\mathrm{dt}}=k_{\mathrm{aTc}}\left(u_{\mathrm{a}}-aT{c}^i\right)$$

$$\frac{\mathrm{d}IPT{G}^i}{\mathrm{dt}}=k_{\mathrm{IPTG}}\left(u_{\mathrm{p}}-IPT{G}^i\right)$$

$${\Phi }_{\mathrm{T}}(t):=\frac{1}{1+{\left(\frac{Tet{R}^i}{{\theta }_{\mathrm{TetR}}}\cdot \frac{1}{1+{\left(\frac{aT{c}^i}{{\theta }_{\mathrm{aTc}}}\right)}^{{\eta }_{\mathrm{aTc}}}}\right)}^{{\eta }_{\mathrm{TetR}}}}$$

$${\Phi }_{\mathrm{L}}(t):=\frac{1}{1+{\left(\frac{Lac{I}^i}{{\theta }_{\mathrm{LacI}}}\cdot \frac{1}{1+{\left(\frac{IPT{G}^i}{{\theta }_{\mathrm{IPTG}}}\right)}^{{\eta }_{\mathrm{IPTG}}}}\right)}^{{\eta }_{\mathrm{LacI}}}}$$

$$\mathrm{d}\mathbf{x}(t)=S\cdot a(\mathbf{x}(t))\cdot\mathrm{dt}+S\cdot\mathrm{diag}\big{(}\sqrt{a(\mathbf{x}(t))}\big{)}\cdot\mathrm{d}\mathbf{w}$$

$$\underset{t\to \infty }{\lim }|e_{\mathrm{A}}(t)|<ϵ\text{ and }\underset{t\to \infty }{\lim }|e_{\mathrm{B}}(t)|<ϵ,$$

$$u_{\mathrm{p}}=\left(1-\frac{u_{\mathrm{a}}}{U_{\mathrm{a}}}\right)U_{\mathrm{p}}$$

$$u(t)=\{\begin{array}{ll}{u}_1,& |e_{\mathrm{B}}(t)|\ge |e_{\mathrm{A}}(t)|\\ u_2,& |e_{\mathrm{B}}(t)|<|e_{\mathrm{A}}(t)|\end{array},$$

$$u_1=\{\begin{array}{ll}{\left[0,U_{\mathrm{p}}\right]}^{\top },& e_{\mathrm{B}}(t)\le 0\\ {\left[U_{\mathrm{a}},0\right]}^{\top },& e_{\mathrm{B}}(t)>0\end{array},u_2=\{\begin{array}{ll}{\left[U_{\mathrm{a}},0\right]}^{\top },& e_{\mathrm{A}}(t)\le 0\\ {\left[0,U_{\mathrm{p}}\right]}^{\top },& e_{\mathrm{A}}(t)>0\end{array}.$$

$$\begin{array}{rl}{u}_{\mathrm{a}}(t)=& k_{\mathrm{P},\mathrm{a}}{e}_{\mathrm{B}}(t)+k_{\mathrm{I},\mathrm{a}}{\int }_0^t{e}_{\mathrm{B}}(t)dt\\ & -\left(k_{\mathrm{P},\mathrm{p}}{e}_{\mathrm{A}}(t)+k_{\mathrm{I},\mathrm{p}}{\int }_0^t{e}_{\mathrm{A}}(t)dt\right)\end{array}$$

$$\{\begin{array}{ll}{u}_{\mathrm{a}}\in \left[0,50\right],& \text{if }|e_{\mathrm{B}}|<|e_{\mathrm{A}}|\\ u_{\mathrm{a}}\in \left[0,100\right],& \text{otherwise}\end{array}$$

$$J(\text{x},r,u,t)={\int }_0^{Tp}\left(\alpha \Vert e_{\mathrm{B}}(t)\Vert +(1-\alpha )\Vert e_{\mathrm{A}}(t)\Vert \right)dt$$

Full text

I Introduction

The aim of Synthetic Biology is to engineer biomolecular systems to achieve new useful functionalities [1]. Potential applications range from designing bacteria that can produce biofuels or sense and degrade pollutant in the environment (like hydrocarbons and plastic), to immune cells that can track and kill cancer cells, or that can release drugs at specific points and conditions to avoid side effects (see [1] for references). This is possible by designing genetic circuits with programmed functionalities and embedding them into living cells. However, most of the engineered genetic circuits have been designed to work at single-cell level. As a consequence, their functional complexity is limited by inherent factors such as excessive metabolic burden on the cell, competition of limited resources, and incompatible chemical reactions.

A promising approach to overcome these issues is to engineer synthetic microbial consortia in which the effort is divided and assigned to different subpopulations of cells to achieve more sophisticated functionalities [2]. Recent cooperative consortia designs include a predator-prey system [3], an emergent oscillator [4], a toggle-switch implemented across two species [5], and a multicellular feedback control scheme where the control functions are split between two species [6]. Unfortunately, the correct functioning of a multicellular consortium requires cocultivating and maintaining multiple cell populations. As different cells in the consortium embed specific sets of genetic circuits, they also present different growth rates due to uneven metabolic burdens and might show additional undesired dynamics, such as oscillations [5]. Therefore, when different strains are mixed together, it is essential to maintain their stable coexistence by controlling their relative population numbers (i.e. their ratio). This is usually achieved by encoding in the synthetic design some dynamic equilibrium between the two populations, e.g. [7],[8]. However, if one of the two populations eventually dies out, these solutions can either lead to uncontrolled growth of one population or to the extinction of both. Moreover, the steady-state value of the populations ratio is hard-coded into the genes without any possibility of being changed online.

In this paper we present an alternative approach to control the populations’ ratio in mono-strain consortia by means of external control inputs. Specifically, we consider the case in which there exists a bistable memory mechanism inside the cells, such as the genetic toggle switch circuit [9], whose current internal state defines which of the two possible roles, or “working-condition”, the cell is playing in the consortium. We assume that, by changing the concentrations of some inducer molecules in the growth medium, it is possible to make cells switch their role and hence keep the populations ratio to a desired value. Albeit requiring a possibly more complex design with respect to other multi-strain scenarios, this approach has the advantage of being intrinsically robust to extinction events that could undermine the operation of the entire consortium. Also, it allows the ratio of the populations to be changed online, in real-time, if needed.

The crucial problem we address in this paper is to design feedback control strategies able to steer the inducer molecules inputs to achieve and maintain a desired cell ratio. We define this problem as “ratiometric” control of cell populations.

We propose and test three different external control strategies to regulate the populations ratio to any value, namely a Bang-Bang controller, a PI controller and a model predictive controller (MPC). The control laws are designed taking into account realistic physical and technological implementation constraints that are present in microfluidic-based experimental platforms. Finally, the proposed controllers are validated in-silico using realistic stochastic agent-based simulations in BSim [10] that appropriately model also spatial and diffusion effects in the microfluidic device and cell growth.

II The ratiometric control problem

We want to design some feedback control strategy such that, by acting on some input signals common to every cells, the ratio between the two populations is asymptotically regulated to some desired value.

II-A Population model

We assume that the bistable memory required by the cells to guarantee their the correct operations in the consortium is realized by means of an inducible genetic toggle switch [9]. This genetic regulatory network consists of two repressor proteins, LacI and TetR, both repressing each other’s promoter, so that only one protein is fully expressed at any time. From a modelling viewpoint, the genetic toggle switch is a bistable dynamical system, possessing two stable equilibria, A and B, each associated to a fully expressed protein, and a saddle equilibrium point, whose stable manifold is the boundary separating the regions of attraction of the other two. Thus, given an initial condition, its solutions will converge to one of the two stable equilibria. The expression level of the two repressing proteins can be flipped by changing the concentration of two inducer molecules, aTc and IPTG. This causes the occurrence of two saddle-node bifurcations yielding the required reversible bistable memory function. We use the inducible toggle switch model described in [11] and further analyzed in [12, 13]. Namely, we assume the dynamics of the $i$-th cell in the consortium can be written as follows:

[TABLE]

where the state variables denote concentrations of molecules inside the cell. The parameters $\kappa_{\mathrm{L/T}}^{\mathrm{m0}}$, $\kappa_{\mathrm{L/T}}^{\mathrm{m}}$, $\kappa_{\mathrm{L/T}}^{\mathrm{p}}$, $\gamma_{\mathrm{L/T}}^{\mathrm{m}}$, $\gamma_{\mathrm{L/T}}^{\mathrm{p}}$, $k_{\mathrm{aTc/IPTG}}$ are leakage transcription, transcription, translation, mRNA degradation and protein degradation rates, and diffusion rates of the inducers across the cell membrane, respectively. The variables $u_{\mathrm{a}}$ and $u_{\mathrm{p}}$ denote the concentrations of the inducer molecules in the growth medium and they also represent the control inputs common to every cell in the populations. Moreover, in the previous equations, the input effects are modelled by the following terms:

[TABLE]

All parameter values of are provided in Table I and are also the same used in [11].

The previous dynamical model is a deterministic description of the evolution of the molecule concentrations in the system and, therefore, it is only an approximation of the stochastic biochemical processes taking place inside the cells. To obtain a more accurate description of the stochastic processes governing the dynamics, for validation we adopted the SDE-based algorithm described in [14] that provide a better approximation of the Chemical Master Equation [15, 16] of the system. Formally, we solved:

[TABLE]

where $\mathbf{x}(t)$ is the state of the process, $S$ is the stoichiometric matrix, $a(\mathbf{x}(t))$ is a vector containing the propensity functions associated to each reaction and $\mathbf{w}$ is a vector of independent standard Wiener processes. Both $S$ and $a(\mathbf{x})$ are the same used in [11].

As shown later in the in-silico validation of the control approaches, the heterogeneity in the response of the cells, provided in our case by the biochemical noise, is a fundamental ingredient to solve the ratiometric control problem.

II-B Problem Statement

We denote by $\mathcal{N}_{t}$ the finite set of all cells in the consortium at time $t$, and with $N(t)=|\mathcal{N}_{t}|$ its cardinality. Note that the number of cells may vary in time as a consequence of cell births and deaths or of their accidental outflow from the microfluidic chamber in which they are hosted. We define the following sets: $\mathcal{A}_{t}:=\{i\in\mathcal{N}_{t}:\,TetR^{i}(t)>2\,LacI^{i}(t)\}$, $\mathcal{B}_{t}:=\{i\in\mathcal{N}_{t}:\,LacI^{i}(t)>2\,TetR^{i}(t)\}$, and $\mathcal{C}_{t}:=\{i\in\mathcal{N}_{t}:\,i\notin\mathcal{A}_{t},i\notin\mathcal{B}_{t}\}$. We also denote with $n_{\mathrm{A}}(t)$ and $n_{\mathrm{B}}(t)$ the cardinalities of $\mathcal{A}_{t}$ and $\mathcal{B}_{t}$ at time $t$, respectively. It is clear from Figure 1 that these sets are disjoint and form a partition of $\mathcal{N}_{t}$ for all $t$.

Noticing the relative position in state space $(LacI^{i},TetR^{i})$ of the stable equilibria A and B of the toggle switch, we say that at time $t$ cell $i$ belongs to the population A (population B), if $i\in\mathcal{A}_{t}$ ($i\in\mathcal{B}_{t}$, respectively). Moreover, we define as $r_{\mathrm{A}}(t)=\frac{n_{\mathrm{A}}(t)}{N(t)}$ and $r_{\mathrm{B}}(t)=\frac{n_{\mathrm{B}}(t)}{N(t)}$ the ratio of cells that belong to population A and population B, respectively.

Definition**.**

Given a consortium of cells whose dynamics is described by (1)-(6) and a desired ratio $r\in[0,1]$ of cells belonging to one population, for example B, we say that the control law $u(t)=\left[u_{\mathrm{a}}(t),\,u_{\mathrm{p}}(t)\right]^{\top}$ solves the ratiometric control problem if, for some small positive constant $\epsilon$,

[TABLE]

*where $e_{\mathrm{B}}(t)=r-r_{\mathrm{B}}(t)$ and $e_{\mathrm{A}}(t)=(1-r)-r_{\mathrm{A}}(t)$. *

III Proposed Control Strategies

In this section we propose three control strategies to solve the ratiometric control problem. Specifically, we present a Bang-Bang controller, a PI controller and a model predictive controller (MPC). All controllers are ad-hoc implementations that take explicitly into account the physical and technological constraints related to an experimental microfluidic platform. Then, in Section IV we validate the controllers in-silico. In contrast to [11] where the objective of the controllers proposed was to regulate at an intermediate level the expression of the genes of a single toggle-switch, here the feedback loop is closed on the entire cell population and the cells are split in two groups, each fully expressing either of the two genes.

Experimental Constraints

The experimental platform we consider as a reference is based on microfluidics as the one described in [17, 18] which uses a fluorescence microscope to measure the current state of the cells. We then have to take into account the following realistic constraints:

the state of the cells cannot be measured more often than $5\,\mathrm{min}$ to avoid excessive phototoxicity; 2. 2.

there is a time delay between $20$ and $40\,\mathrm{s}$ on the actuation of the control inputs due to the time that the flow of the chemical inducers takes to reach the chambers on the microfluidic chip where cells are hosted; 3. 3.

the minimum time interval between two consecutive control inputs cannot be less than $15\,\mathrm{min}$ to limit excessive osmotic stress on the cells; 4. 4.

the maximum duration of any experiment cannot exceed $24$ hours, to avoid substantial cell mutations during the experiments.

Moreover, the specific implementation of the microfluidic device also introduces constraints on the possible classes of input signals $u(t)=\left[u_{\mathrm{a}}(t),\,u_{\mathrm{p}}(t)\right]^{\top}$ that can be generated by the actuators.

We consider two possible implementations:

a T-Junction, which limits $u_{\mathrm{a}}$ and $u_{\mathrm{p}}$ to be mutually exclusive and with fixed amplitudes, that is $u$ is either set to $[U_{\mathrm{a}},0]^{\top}$, which causes $e_{\mathrm{B}}$ to decrease and $e_{\mathrm{A}}$ to increase, or to $[0,U_{\mathrm{p}}]^{\top}$, which does the opposite; 2. 2.

a Dial-A-Wave (DAW) system [19], which constraints $u_{\mathrm{a}}$ and $u_{\mathrm{p}}$ to be in a convex combination. Namely, given $u_{\mathrm{a}}\in[0,U_{a}]$ we have

[TABLE]

In the above equations, $U_{\mathrm{a}}\in\left[0,100\right]$ and $U_{\mathrm{p}}\in\left[0,1\right]$ are control amplitudes to be selected and denote the maximum concentrations possible of the inducers that are present in the reservoirs (These values are the same as those that were used in-vivo in [11].).

Depending on which implementation is considered, only specific controllers are feasible. More precisely, for the T-Junction implementation only a Bang-Bang controller is considered, while for the Dial-A-Wave system we design a PI controller and an MPC.

III-A Bang-Bang Controller

The Bang-Bang controller implemented via a T-junction consists of two mutually exclusive inputs with fixed amplitude which are applied to the system depending on the current value of the error signals $e_{\mathrm{A}}(t)$ and $e_{\mathrm{B}}(t)$. Specifically, at any time $t$ the input is applied that causes the $\max\{|e_{\mathrm{A}}(t)|,|e_{\mathrm{B}}(t)|\}$ to decrease. More formally, the control input $u(t)=\left[u_{\mathrm{a}}(t),\,u_{\mathrm{p}}(t)\right]^{\top}$ is chosen as

[TABLE]

where

[TABLE]

III-B PI Controller

The PI control inputs to be implemented via the Dial-a-Wave system are chosen as:

[TABLE]

with $k_{\mathrm{P,p}}$, $k_{\mathrm{P,a}}$, $k_{\mathrm{I,a}}$ and $k_{\mathrm{I,p}}$ being the control gains, and $u_{\mathrm{p}}(t)$ given by (9).

Moreover, to improve the performance and guarantee that the control signals do not exceed their admissible values, the PI controller is complemented by a dynamic saturation defined as:

[TABLE]

and an anti wind-up scheme.

III-C MPC Algorithm

The last algorithm we considered is a Model Predictive Controller (MPC) [20]. Given the state of the cell population, say $\textbf{x}=[x_{1},\dots,x_{N(t)}]^{\top}\in\mathbb{R}^{2N(t)}$, where $x_{i}=[LacI^{i},TetR^{i}]^{\top}$ is the state of the $i$-th cell, we compute the optimal control input over the time interval $\left[t,t+T_{p}\right]$ which minimizes the cost function

[TABLE]

with $T_{p}$ being the controller prediction time and $\alpha\in(0,1)$ a constant design parameter. We then apply the computed optimal control input over the interval $\left[t,t+T_{c}\right)$ where $T_{c}<T_{p}$ is a control time to be selected during the implementation.

To reduce the computational burden, in our implementation we evaluated the cost function $J(\textbf{x},u)$ using only a representative subset of cells chosen as a sample of the entire population. This subset of cells is chosen such that the reduced ratios (i.e. the ratios $r_{\mathrm{A}}$ and $r_{\mathrm{B}}$ computed on the subset) is as close as possible to the ratios of the entire population. Also, a genetic algorithm taken from [21] (without the mutation phase) was used to find the optimal control sequence.

IV In-silico validation and comparisons

We tested all the designed control laws in-silico. First, we performed a batch of simulations in Matlab assuming a constant population of 30 cells, then we performed more accurate simulations using an agent-based simulator specifically designed for bacterial populations called BSim [10].

In all simulations we consider a desired ratio $r=0.6$, with initial conditions taken in a neighborhood of the saddle point between A and B, specifically we picked the initial conditions for each cell with a random uniform distribution in the intervals: ${mRNA_{\mathrm{LacI,0}}^{i}}\in\left[3,6\right]$, ${mRNA_{\mathrm{TetR,0}}^{i}}\in\left[3,6\right]$, ${LacI}_{0}^{i}\in\left[150,300\right]$, ${TetR}_{0}^{i}\in\left[200,400\right]$. To mimic the real experimental constraints, the state of the population is sampled every $T_{s}=5\,\mathrm{min}$ and the control inputs are evaluated every $T_{c}=15\,\mathrm{min}$.

IV-A Numerical simulations in Matlab

For the Bang-Bang controller we empirically set the control amplitudes to $U_{\mathrm{a}}=60$ and $U_{\mathrm{p}}=0.5$. With these values we obtained the evolution of errors and control inputs shown in panels (a) and (b) of Figure 2. We notice that the Bang-Bang controller achieves good performance with a settling time of about 1300 min.

For the PI controller, the control input amplitudes and control gains were empirically selected as $U_{\mathrm{a}}=100$, $U_{\mathrm{p}}=1$, $k_{\mathrm{P,a}}=66.67$, $k_{\mathrm{P,p}}=2.25$, $k_{\mathrm{I,a}}=1.2$ and $k_{\mathrm{I,p}}=0.006$ obtaining the results portrayed in panels (c), (d) of Figure 2. In this case we observe a settling time which is half the one obtained with the Bang Bang controller, together with lower error values at steady state.

Finally, for the MPC controller, we chose $[U_{\mathrm{a}},U_{\mathrm{p}}]=[60,0.5]$, $T_{p}=75\,\mathrm{min}$, $T_{c}=T_{s}=15\,\mathrm{min}$, $\alpha=0.6$. The genetic algorithm’s parameters were set to $N_{\mathrm{p}}=20$, $M_{\max}=10$, where $N_{p}$ is the length of the control sequence generated at each step and $M_{\max}$ the number of generations being considered.

The resulting errors and the control inputs are shown in panels (e) and (f) of Figure 2, which confirm the MPC as the best strategy with a settling time almost $40\%$ shorter than the one observed with the PI.

IV-B Agent-based simulations in BSim

To provide a more realistic validation of our control strategies, we used BSim [10], a realistic agent-based simulator of bacterial populations, which considers also cells reproduction, spatial distribution and geometry, the diffusion of the chemicals in the environment and, more importantly, the flush-out of cells from the chamber. To run the required stochastic simulations, we extended BSim with an Euler-Maruyama solver [22].

As a reference for the geometry of the microfluidic device we used a scaled version of the one described in [23] with dimensions $13.3\mu\mathrm{m}\times 16.6\mu\mathrm{m}\times 1\mu\mathrm{m}$ which can contain a population of about $50$ cells.

We tested in BSim all the proposed controllers and snapshots of a typical BSim simulation is shown in Figure 3 where cells fully expressing one of the two repressor genes are depicted either in red or green.

The errors obtained via BSim in-silico experiments are reported in Figure 4. It can be noticed that the fluctuations of the error signals are higher than in the previous Matlab simulations essentially due to cell growth and splitting, and the flush-out of the cells from the chamber. However, the average error evolution is qualitatively the same, confirming the good performance of the controllers.

IV-C Performance Comparison

To compare the performance of the controllers, we considered the following performance indices evaluated averaging over $M$ simulation trials:

(i) the average value of the error norm over the total simulation time $T_{\mathrm{sim}}$, $\bar{e}=\frac{1}{M}\sum_{j=1}^{M}\left(\frac{1}{T_{\mathrm{sim}}}\int_{0}^{T_{\mathrm{sim}}}\lVert e_{j}(t)\rVert dt\right),$ and (ii) over the last $180\,\mathrm{min}$, $\bar{e}_{\mathrm{f}}=\frac{1}{M}\sum_{j=1}^{M}\left(\frac{1}{180}\int_{T_{\mathrm{sim}}-180}^{T_{\mathrm{sim}}}\lVert e_{j}(t)\rVert dt\right),$ where for the $j$-th trial $e_{j}(t)=[e_{\mathrm{A}}^{j}(t),\,e_{\mathrm{B}}^{j}(t)]^{\top}$,

and (iii) the average settling time at $15\%$ of the error, $\bar{t}_{\mathrm{s}}$.

In Table II we report the values of the previous indices considering $M=30$ for the simulations in Matlab and $M=1$ for those in BSim. It can be observed that the MPC algorithm guarantees excellent performance both in terms of settling time and steady-state error norm. Therefore, it is the best candidate for the in-vivo implementation that will be the next stage of our ongoing research. Nevertheless, the Bang-Bang controller offers a good compromise between ease of implementation and performance.

V Conclusions

We considered the ratiometric control problem in a mono-strain microbial consortium made of bacteria embedding a bistable toggle switch. We demonstrated that, by varying global inputs to all the cells, it is possible to control the ratio between those stabilizing onto one equilibrium and those on the other. Namely, we presented three control strategies to regulate the cells in the consortium to the desired ratio. The control design took into account the constraints of a possible experimental microfluidic implementation. We tested performance of the controllers in-silico by numerical and realistic agent-based simulations. In both cases, it emerged that the MPC algorithm guarantees excellent performance in terms of both the settling time and the steady-state error values. Future work will be aimed at validating in-vivo the proposed controllers and exploiting them for multicellular feedback control schemes such as the one described in [6].

ACKNOWLEDGMENT

The authors wish to acknowledge support from the research project COSY-BIO (Control Engineering of Biological Systems for Reliable Synthetic Biology Applications) funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 766840.