Modelling Interactions of a Microbial Community through an Optimal Tracking Control Approach

Genetic sequencing measurements of the existing species in a reactor are incorporated into a classic chemostat model through a non linear optimal tracking approach. The framework is applied to data coming from a nitrification process, resulting in a dynamical system with a high dimensional state. The observed ecological dynamics from the bacterial species are able to explain the substrates dynamics.


I. INTRODUCTION
In the last decades genetic sequencing techniques have become an affordable measurement in bioprocesses.However using this measurements for the prediction and control of the bioreactors performance (in terms of proportions of chemicals consumed and produced) lacks a modelling framework [1].
It has been theorised that interactions between microorganisms could help explain the difference in performance of reactors operating under very similar environmental conditions [2].Some types of interactions have already been accounted in chemostat theory, such as inhibition in the presence of other species, or density dependence.Nevertheless this represent just a fraction of the possible ecological interactions [3] that could take place in a microbial community.
The central idea behind the method presented in this work is that if one can accurately follow the dynamic of the species present in an ecosystem for a certain time horizon then one can recover the underlying interactions that might be driving the system.

II. MODEL DEFINITION AND PROBLEM STATEMENT
The following conventions are used: (1) Let m be an integer then [m] := {1, . . ., n}, (2) For integers m 1 and m 2 and a ∈ R, a m 1 ×m 2 represents a matrix of m 1 rows and m 2 columns with a in every entry, (3) let m be an integer then I m is the identity matrix of size m, (4) given a vector v = (v 1 , . . ., v n ) ∈ R n the function diag(v) returns a square matrix M with entries M ii = v i , ∀i ∈ [n], and M i j = 0, ∀i, j ∈ [n], i = j, and (5) Let a 1 , a 2 ∈ R then a 1 ∨ a 2 := max{a 1 , a 2 } and a 1 ∧ a 2 := min{a 1 , a 2 }.

A. Stoichiometric Equations
Consider the situation where n different Operative Taxonomic Units (OTU) are present in a system.A group of OTU (G 1 ⊂ [n]) consumes a substrate s 1 for growth and as a result s 2 is produced, while another group of OTU (G 2 ⊂ [n]) consumes s 2 for growth and produces s 3 .G 1 and G 2 are called functional groups.This situation is often referred to 1 LBE, INRAE, Univ.Montpellier, 11100, Narbonne, France.* Corresponding author : pablo.ugalde-salas@inrae.fr as a cascade bioreaction.The number of organisms in G 1 and G 2 will be denoted n 1 and n 2 , respectively.It will be assumed that sets G 1 and G 2 satisfy: The situation is schematically described in bioreactions (R G1) and (R G2).The terms y i are known as yields, they represent the number of moles of biomass produced per mole of substrate consumed.The term M x represents a molecule of biomass and several expressions can be found in the literature (e.g.CH 1.613 O 0.557 N 0.158 [4]).s = (s 1 , s 2 , s 3 ) and x = (x 1 , . . ., x n ) denote the vectors of concentrations of substrates and of OTU present in the system, respectively.Each reaction has its own process rate (also known as growth function or kinetics) µ i (s, x).Note that it depends on x to show the influence of other OTU on the growth function.A classic treatment of general reactions and kinetics in biotechnology can be found in [5].

B. Mass Balance Equations
The system where reactions (R G1) and (R G2) take place is a container of fixed volume V containing a liquid medium with the same input and output flow Q = Q(t), which can vary in time.The dilution rate is then defined as D = D(t) = Q V .The input flow contains a time dependent concentration s in = s in (t) of the same chemical composition than s 1 , for readability the time dependence notation is dropped.For chemostat theory detailed modelling and mathematical analysis one can refer to [6].The dynamics of the described situation are seen in (1).

C. Kinetics and Interactions
Given a vector (v 1 , . . ., v n ) the interaction function I is noted as: The growth function of each OTU will be assumed to be a product of substrate dependent function and the interaction function.Let f i (s) be a bounded, positive, and continuous function of s (e.g.Monod, Haldane).The growth equation of OTU i becomes: The notation f (s) := ( f 1 (s), . . ., f n (s)) will be used.
Since the growth of a single strain in chemostat experiments is limited by the substrate concentration , when no interactions are present one should recover expression f i (s).Therefore if there are no interactions then I i (x) = 1.From this reasoning, note that lim x→0 I i (x) = 1 since if there is minimal presence of OTU, interactions can not exist.Furthermore it is assumed that I i (•) is a continuously differentiable function on x.For making explicit all of the former: The interaction function I previously defined satisfies: By forcing the positivity of the interaction function then the growth functions are still positive and the system remains well posed [6], meaning that if system (1) is provided with non-negative initial conditions, the trajectories remain nonnegative.
The problem is then to reconstruct the interaction function from the observations of the microbial community composition in time.

D. Optimal Control Setting for Unravelling the Interaction Function
Suppose that the functions f i (s), and the yields y i are well-known.By using experimental measurements of x, represented by z(t), the objective is to reconstruct function I(x).For doing so, the terms I i (x) are replaced by a control u i (t) (u i (t) = I i (x(t))) and a state feedback control law is obtained by solving a nonlinear optimal tracking problem.Consider the observable system (4).
Let Q, and R be positive definite matrices, and let • Q and • R be the weighted norms by the corresponding matrices.The optimal tracking problem is defined as: The control u(t) is intended to drive the system to be near a desired output z(t), which in this context are the measurements of the concentrations of OTU.The term (u − 1 n×1 ) R , was added for two reasons: First because the interest is testing the idea that interactions could be driving the system.Therefore adding a penalization in the objective function for each control to remain near 1 can be seen as an attempt to explain data without any interaction.
Second, to force a regularized control.Otherwise note that u is linear in the state equations (4), therefore if the integral cost does not have a non-linear expression of u the optimal control will be of a bang-bang type with possibly singular arcs.Since the objective is to find a differentiable expression of I(x) the addition of the regularization term is deemed necessary.
The problem of approximating the solution of the system to a desired reference (z in this case) is called the optimal tracking problem.For solving such a problem the approach developed by Cimen et al. [7] was adapted to our problem, using [8].The method proposed involves the resolution of Approximating Sequences of Ricatti Equations.
The method consists on iteratively calculating trajectories of system (4) with a certain control and feeding a nonautonomous Ricatti differential equation with the resulting trajectory.Then a control law that uses the Ricatti equation is proposed and a new trajectory is calculated.The iteration is stopped when a convergence in the output or in the control is observed.
First the tracking problem does not consider a constrained control, so one should consider a reformulation as proposed by the same author in other works [8].Nevertheless, the methods in [7] were directly used, and an explicit constraint in the synthesis of the control was added for the control to remain within the desired bounds (see equation ( 21)).
Even though this is probably suboptimal it worked for our purposes.
Second in their problem they used an approximation of their original dynamics for calculating the trajectories and subsequently proved that such linearisation converged to the original dynamics.In the case here presented the linearisation of the dynamics was not necessary.

E. Tracking Problem Reformulation
For applying the methods developed in [7].Define the system state X = (x, s), and the change of variables v i = u i −1 with v = (v 1 , . . ., v n ).System (4) may be rewritten then as: Define A 12 (X) = 0 n×3 (13) Then system (4) becomes: Let z(t) ∈ R n be the vector containing the experimental measurements of the OTU concentrations in time.The cost functional is now rewritten as (19).
Define the dynamic sequences for i ∈ N, X [i] as : (20) And for i = 0 let X [0] (t) = X 0 .The control law is given by Where P [i] (t) ∈ M n+3×n+3 (R) and s [i] f (t) ∈ R n+3 are the solution to the systems of differential equations ( 22) and (23), respectively.
Note that system ( 22) and (23) consists of (n + 3) 2 + (n + 3) differential equations.In the following the inspection of system ( 22) and ( 23) allows us to reduce the number of differential equations to 2n.Replacing matrices ( 9), (10), and (17) in ( 22) and ( 23) one obtains: The fact that certain entries are autonomous in the dynamic and that the constantly zero function is a solution for them, implies by existence and uniqueness that they should be constantly zero.Then P [i] has n × n non zero entries and s [i] has n non zero entries, explicitly: Where P[i] (t) and s f [i] (t) satisfy differential equations (28) and (29), respectively.Note that matrices depend on the previous state iteration i.e.A 11 = A 11 X [i−1] (t) and Assume now that matrices Q and R are diagonal matrices.Inspecting the equation (28) one notices that if i = j, P i j (t) = 0 is a solution of all non diagonal entries.And therefore, once again, by existence and uniqueness they are constantly zero.Hence in system (28) only the diagonal entries should be calculated as seen in equation (30) .
And the control law is given by Differential equations ( 30) and (29) were solved using standard backward numerical integration for each iteration.Note that from the former computations the method relies in calculating 2n differential equations instead of (n+3) 2 +(n+ 3) as the straightforward application of the method would have meant.

III. APPLICATION AND DISCUSSION
The choice of matrices Q and R are always of discussion in a quadratic regulator.Q represents the covariance matrix of measurements, while the interpretation of R is not clear.In the simulations here presented Q = I n and R = 0.0001I n , found by trial and error in order for the control to fit accurately the species abundance.The method stopped after 110 iterations where each iteration took approximately 600 seconds of computing time in a computer equipped with 8gb of RAM memory and Intel core i3-7100U CPU 2,40 GHz.
The tracking problem was applied to data coming from a nitrification process; the experimental conditions are described in [9].For exploring the hypothesis of interactions as drivers of bioreactors performance environmental conditions should be kept as constant as possible.Therefore only data from day 183 onwards was used because a change in the operating temperature happened at that point, which is known to have an effect on the kinetics.For choosing which species belong in which functional group, the procedure described in [10] was used: from day 183 to day 230: 17 OTU were identified in the G 1 group (known as Ammonium oxidizing bacteria AOB) and 8 in the G 2 functional group (known as Nitrite oxidizing bacteria NOB).The system state dimension becomes 28.
The knowledge of functions f i (s) was based on a study of nitrification's kinetic parameters [11].Particularly given the system's ammonium and nitrate concentration a Monod function (eq (32)) was used for G 1 and G 2 .
For the reader to gain understanding of the situation, a simulation of the system using the experiments operating parameters (D and s in ) is presented without control (i.e.v(t) = 0) in Figure 1.s 3 accumulates all along the trajectory, but when compared to data it is clear that s 3 stops accumulating after a while.
When applying the presented method one obtains the simulation that can be seen in figure 2. The method captures the tendencies in the measured substrates, but not their absolute values.The simulation of each functional G 1 (AOB), and G 2 (NOB) can be seen in figures 3 and 4, respectively.Note that the fit of G 1 seems better than the fit of G 2 , which can be explained by the 1 order of magnitude difference of abundance between the members of the group.One can observe that in figure 2 the predicted s 3 around day 240 is much higher than the measurement.If one looks at figure 4 it suggests attributing this extra production of s 3 around day 240 to the growth of OTU 24.What would have happened if OTU 24 had been classified in the other functional group?One would have seen a decrease in s 1 , an increase in s 2 , and a decrease in s 3 which could have better fitted the measurements.Thus suggesting that OTU 24 was misclassified from the method presented in [10].
In figures 5 and 6 it can be seen the obtained controls for each functional group.Both controls saturate at certain points, but in general terms they vary throughout the dynamic.
Recalling that u(t) is supposed to represent I(x(t)) (the interactions terms), one could try to fit different forms of I(•) for example I(x) = 1 + Ax, where A is a real valued square matrix.In figure 7 the simulation of the system after the identification of matrix A is presented, suggesting that an affine interaction function is not able to capture the behaviour of the system.

IV. CONCLUSION AND PERSPECTIVES
Considering the dimension of the system analysed, the approach here presented seems satisfactory for the integration of genetic sequencing in the context of bioreactors, particularly the Ricatti equations grow linearly in the number of species, rendering the method scalable for a bigger number of species.However some issues should be taken care of: 1) test the method on synthetic data for a class of interaction functions 2) the optimality of the control should be improved, considering is just a truncated control from a state feedback 3) A clear rational on the choice of matrices Q and R. Finally the hypothesis of misclassification presented for OTU 24 proposes as well a retroactive way to reclassify an OTU if the substrates predictions go astray.

Fig. 7 :
Fig. 7: Simulation obtained for the interaction function I(x) = 1 + Ax, where A was fitted with the obtained controls from figures 5 and 6.