Control of multidimensional systems on complex network

Multidimensional systems coupled via complex networks are widespread in nature and thus frequently invoked for a large plethora of interesting applications. From ecology to physics, individual entities in mutual interactions are grouped in families, homogeneous in kind. These latter interact selectively, through a sequence of self-consistently regulated steps, whose deeply rooted architecture is stored in the assigned matrix of connections. The asymptotic equilibrium eventually attained by the system, and its associated stability, can be assessed by employing standard nonlinear dynamics tools. For many practical applications, it is however important to externally drive the system towards a desired equilibrium, which is resilient, hence stable, to external perturbations. To this end we here consider a system made up of N interacting populations which evolve according to general rate equations, bearing attributes of universality. One species is added to the pool of interacting families and used as a dynamical controller to induce novel stable equilibria. Use can be made of the root locus method to shape the needed control, in terms of intrinsic reactivity and adopted protocol of injection. The proposed method is tested on both synthetic and real data, thus enabling to demonstrate its robustness and versatility.


I. GENETIC NETWORK MODEL
We shall here justify the model of genetic regulatory network analyzed in the main text. Consider first a small regulatory network consisting of one gene (whose activity is labelled x) and one protein (associated to the continuous concentration y). A positive regulation loop can be modeled as: x = k 1 g(y) − γ 1 ẋ y = k 2 x − γ 2 y where: g(y) = y n K + y n . (S1) In the following, we will make the choice n = 2 and K = 1. Similarly, a negative regulation loop can be modeled as:ẋ = k 1 (1 − g(y)) − γ 1 ẋ y = k 2 x − γ 2 y As it should be, the concentration of proteins grows with the level of gene activity. As a first step approximation, and to avoid dealing with two distinct families of mutually interlinked constituents (genes and proteins), we can replace y with x in the argument of g(·), via adiabatic elimination (apart from redefinition of the involved constants). Building on the above we model the extended regulatory network as: where x i stands for the activity of gene i. For positive feedbacks between species i and j, A ij = 1, while, for negative loops A ij = −1. The parameter η i stands for the number of negative loops (number of negative entries of the i-th row of A) that are associated to node i. Finally, to keep the structure as simple as possible we set k i = 1 for all i. The matrix A employed in the example reported in the main body has a simple structure and it has been chosen for purely illustrative purposes: it represents a regular tree network with branching ratio r = 4.

II. CONTROLLING THE MICROBIOTA NETWORK: THE EXPERIMENTAL DATA
The second application of the proposed control method deals with a model of gut microbiota. As explained in the main body of the paper, the model is well established and builds on a generalization of the celebrated Lotka-Volterra equations:ẋ The above equations take into account species-species interactions and self-regulation, this latter effect being described in terms of a logistic growth. Despite its intrinsic simplicity, the model is often invoked to explain how ecological interactions, e.g., mutualism and competition for nutrients, can lead to complex phenomena, as multi-stability or antibiotic mediated catastrophic losses of biodiversity. In [? ] the model has been applied to a relatively small (mice gut) microbiota system made up of 11 distinct populations. More precisely, the ten most abundant species have been identified: all together they account for the vast majority (∼ 90%) of the total populations found in the mice gut. The remaining populations are grouped into a unique (non-homogeneous) category referred to as "Other". The authors of [? ] analyzed the experiments reported in [? ] and provided a quantitative characterization of the coefficients that enter the definition of the relevant quantities r, s and A. These latter are reported in table ( The second application discussed in the main body of the paper, i.e. that targeted to controlling the microbiota dynamics, assumes the above parameters. Notice that the set of considered species includes the spore-forming pathogen Clostridium difficile. To lower its concentration (and so diminish the probability of infection) is one of the goals of the implemented control. As explained in the paper, three different control schemes are considered. In the following we will provide some additional information for each of the analyzed schemes. The physical dimension of the inserted controller u is again 10 11 rRNAcopies/cm 3 . The parameters β have dimension of the inverse of a time, while ρ is a-dimensional.

A. Case A: Stabilizing an unstable fixed point by means of an external controller
The first of the three different control strategies for the microbiota network discussed in the main body involves only 5 out of the 11 populations being examined (more precisely Barnesiella, Blautia, und. Mollicutes, Coprobacillus and und. Enterobacteriaceae). An unstable fixed point of equation (S3), which contains only these latter species, has been calculated and denoted byx (see table (S5)). The maximum eigenvalue of the Jacobian matrix evaluated atx is (λ Re ) max = 0.0148 > 0, thus implying instability. As discussed in the main body of the paper, it is possible to stabilize a different equilibrium, denoted by x * in (S5), which corresponds to a slight modification of the unstable solutionx. Following the procedure detailed in the paper one can readily calculate the parameters α and β (S5). The values obtained are reported in table (S5). Direct simulations of the controlled system, as displayed in Figure  S1, confirms that stability has been indeed achieved.  Figure S1: Numerical integration of the controlled system (see equations (2) in the main body). The equilibrium state stabilized upon injection of the controller (stars) is a slight modified version of the initially unstable fixed point, see table (S5). The system is initialized out of equilibrium and, after a transient, converges to x * .

B. Case B: Acting with one species of the pool to damp the concentration of the pathogens
Consider the stable fixed point x * (table (S6)) composed by Barnesiella, und. Lachnospiraceae, Other, Blautia and C. difficile. As explained in the paper, we now insert in the system another population, selected from the extended pool of interacting species. This latter configures as the controller. Vector α therefore follows from the interaction matrix A. One can then calculate the fixed point that can be eventually attained by the system, given the specific selected controller. Retaining only the meaningful cases (fixed points with all positive entries), we obtain three possible solutions: x * L where the added species is uncl. Lachnospiraceae, x * M adding as external control the und. uncl. Mollicutes and x * E adding und. Enterobacteriaceae (see table S6). From inspection of the obtained solutions, one can appreciate the impact of the different controllers employed: in the latter case the concentration of C. difficile stays almost constant, in the second example it increases, while in the first case it is reduced by a significant amount.

Populations
x

C. Case C: Driving to extinction one species, the other being the target of the control
In this case we aim at enforcing the extinction of one of the species (here C. difficile), acting on the other ones. In other words the parameter vector α (which characterizes the action of the control against the species) is forced to have the component relative to C. difficile equal to zero. At the same time, the corresponding entry of the vector x * to be eventually stabilized is set to a negligible value. The concentrations of the species for respectively the initial and the final fixed points are compared in Figure 2 of the main body. The corresponding values of α are also plotted. As a supplementary material, we here report in Figure S3 the root locus diagram obtained for the case at hand. For a description of the symbols, refer to the caption of Figure S2.

III. EXPLOITING A TRANSIENT CONTROL TO DRIVE THE SYSTEM TOWARDS AN EXISTING STABLE FIXED POINT
As anticipated in the main body of the paper, the control can be effectively employed to steer the system towards a stable fixed point of the unperturbed dynamics, starting from out-of-equilibrium initial conditions. In this case, we setu = −γu − j β j (x j − x * j ). Here, x * j is an equilibrium solution of the uncontrolled dynamics, which proves linearly stable to external perturbations. The parameter γ can be tuned as desired so as to help the convergence towards x * j without falling in the basin of attraction of other existing fixed points. As a proof of principle of the method, we choose a stable fixed point of the global macrobiota ecosystem, i.e. including the complete pool of 11 populations. This is characterized by x * 1 = 9.299, x * 3 = 12.3085, x * 4 = 3.1627 and x * j = 0 for j = 1, 3, 4. The largest real part of the eigenvalues of the associated 11 × 11 Jacobian matrix turns out to be (λ Re ) = −0.1306 < 0, thus implying stability of the aforementioned equilibrium. Imagine to initialize the system out of equilibrium with all species, including the pathogen C. difficile, being assigned a random concentration x j (0) = 0. The system is let evolve for a while and then, at time t * , the control is injected. Here, α and β are assigned as random, uniformly distributed over [0,1], parameters. As clearly displayed in Figure S4, the system is steadily moved towards the equilibrium x * j (stars), while the control converges to zero after an abrupt jump. In other words, after a transient whose duration depends on the chosen parameters, the system achieves its asymptotic (pathogen free) equilibrium and the control can be safely disconnected.

IV. ON THE CONTROLLABILITY CONDITION
Observe that the control procedure here discussed requires computing the vector β. Determining this latter implies inverting a matrix, an operation that imposes a mathematical constraint that we shall hereafter analyze more in depth. The matrix to be inverted, as defined in the main text, reads:  Figure S4: Driving the system towards a stable fixed point of the unperturbed dynamics. The system is initiated out-of-equilibrium: the concentration of all species, including the pathogen C. difficile, is set to values different from zero. At t * = 5 days the control u is injected. After a transient the system converges to its stable equilibrium characterized by x * 1 = 9.299, x * 3 = 12.3085, x * 4 = 3.1627 and x * j = 0 for j = 1, 3, 4, while the control u is turned to zero. Here γ = 10 and u(0) = 15.
where c stands for the coefficients 1 of the characteristic polynomial of G, namely det(G − λI) = N l=0 c l+1 λ l . The invertibility is ensured if det(H) = 0, in formulae: where i1,...i N is the Levi-Civita tensor. This complicated expression can be heavily simplified. The Levi-Civita symbol is in fact totally antisymmetric in the permutation of its indices. As a consequence, all the terms multiplied by i1,...i N which are symmetric in the permutations, cancel out. It follows that the terms containing the product of two or more factors (G k q) with the same power k are identically equal to zero. The only terms which surviveare those obtained by just retaining the largest possible value ofk in each summation (k 1 = N − 1, k 2 = N − 2,..., k N −1 = 1). In formulae: where use has been made of the fact that c N +1 = (−1) N (see equation (S8) in the footnote 1).
Drawing on this preliminary observations it is possible to re-interpret the above controllability constraint, making contact with standard control theory. The controllability condition amounts to require that the matrix has maximum rank. Here, with the notation (G l q) k we identify the k-th entry of the vector obtained from the product of the matrix G to the power of l with vector q. In system theory the matrix C is called the controllability matrix of the pair (G, q). Since C is a square matrix, the maximum rank condition is equivalent to require det(C) = 0, in where with G (j 1 j 1 )(j 2 j 2 )...(j k−1 j k−1 ) we identify the minor obtained from matrix G by removing the j 1 -th, j 2 -th,...,j k−1 -th rows and columns. formulae: where use has been made of the definition of C (S11), namely C ij = (G j−1 q) i . Expressions (S10) and (S12) are identical (except for the sign) and consequently the two conditions, det(H) = 0 e det(C) = 0, prove equivalent. Stated differently, the condition of invertibility of matrix H, obtained as a selfconsistent constraint for the introduced control scheme, coincides with the standard controllability condition, as known in control theory. As a final remark we recall (see main body) that our goal is not to arbitrarily assign the polynomial N (λ) but rather to locate its roots z k within the open left-hand plane. In this respect, the necessary and sufficient systemtheoretic condition is the so-called stabilizability of the pair (G, q), which results in the Popov-Belevitch-Hautus rank condition [? ? ].