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.

Introduction

Investigating the interlinked dynamics of an ensemble composed of units organized in homologous families, constitutes a universal challenge in science, of broad and cross-disciplinary breath [1-3]. Each population is customarily identified in terms of its continuous density. This latter evolves in time, as dictated by specific self-reaction stimuli, that generally bear nonlinear contributions. In a complex and dynamical environment, species experience a large plethora of mutual interactions, declinated via different modalities, notably pairwise exchanges. Cooperative and competitive interference are simultaneously at play, and shape the ultimate fate of the system as a whole [4]. These fundamental ingredients, flexibly combined and properly integrated, are at the roots of any plausible mathematical model targeted to community interactions [5], from ecology [6] to neuroscience [7, 8], passing from genetic and human health [9], through a full load of man-made technological applications [10]. Irrespectively of the specific realm of investigation, each population can be abstractly assigned to a given node of a virtual graph. Directed or indirected edges among nodes exemplify the topological structure of the existing network of interactions [11, 12]. The intricate web of inter-species connections, key information to anticipate the expected dynamics of the system, is therefore encoded in the associated adjacency matrix [13, 14].

In many cases of interest, it is important to drive the system towards a desired equilibrium, that is stable, and thus resilient, to external perturbations [15-19]. For example, hostile pathogens could be forced to go extinct: the stability of the attained equilibrium would efficaciously shield from subsequent harmful invasion and outbreaks. Alternatively, it could prove vital to robustly enhance the expression of species identified as beneficial for the system at hand. Building on these premises, we here develop and test a general control strategy [20-24] targeted to multidimensional systems consisting of a large number of components that interact through a complex network. By inserting one additional species, the controller, which configures as a further node of the collection, we will be able to manipulate the asymptotic dynamics of the system, in terms of existence and stability of the allowed fixed points.

To set the reference frame we will hereafter consider a system consisting of N species (nodes) whose activities = (x1, x2, ⋯, x) obey the coupled nonlinear Eqs [6, 15, 16]:

The first term on the right-hand side specifies the self-dynamics of species i while the second term stems from the interactions of species i with the other species. The nonlinear functions f(x) and g(x, x) encode the dynamical laws that govern the system’s components, while the weighted connectivity matrix captures the interactions between nodes. The elements A can be positive or negative, depending on the specific nature of the interaction, i.e. cooperative or competitive. Notice that system (1) is assumed in [15] as a reference model to analyze resilience patterns in complex networks. Differently from [15], A can here take positive and negative values (see also [25]).

In ecological applications, the number of nodes reflects the biodiversity of the scrutinized sample [6, 16, 26]. Distinct trophic layers materialize as coherent blocks in the adjacency matrix, whose entries modulate the strength of mutual interactions [4]. These are often epitomized by a quadratic response function g(x, x) [16]. Each species is then subjected to a reaction drive f(x), usually a logistic growth with a prescribed carrying capacity [15]. Animals displaying competitive predator-prey interactions or, alternatively, subjected to a symbiotic dependence, such as in plant-pollinator relationships, are among the systems that fall within the aforementioned scenario [4]. Furthermore, the complex community of micro-organisms that live in the digestive tracts of humans and other animals, including insects, can be rooted on similar descriptive grounds [6]. For genetic regulatory networks, the dynamical variables x represent the level of activity of a gene or the concentration of the associated proteins [27]. Species specific reaction terms f(x) account for e.g. degradation or dimerization. The pattern of activation could be effectively modeled by sigmoidal Hill-like functions [1], as follows the classical Michaelis-Menten scheme [28], which incorporates the known map of gene interactions. On a more general perspective, understanding the emerging dynamics in social communities [29], grasping the essence of the learning organization in the brain [30], and implementing efficient protocols for robot navigation in networked swarms [31] are among the very many applications that can be traced back to one of the variants of Eq (1), with a suitable choice of the nonlinear functions f(x) and g(x, x).

Materials and methods

Adding a species to enforce stable equilibria in a multidimensional system

Starting from the above illustrated setting, we will here discuss a suitable control scheme to drive system (1) towards a desired equilibrium , which is linearly stable to externally imposed perturbations. To reach this goal we shall introduce one additional species, the (N + 1)-th component of the collection, suitably designed to yield the sought effect. To set the notation, we indicate by u the component (e.g., concentration, activation level) assigned to the controller and write:

The controller u can exert a direct influence on every component x, as specified by newly added terms α h(x, u) that modify the original system (1). = (α1, α2, …, α) is a vector of N constant parameters, to be self-consistently adjusted following the scheme depicted below. h(x, u) is a generic, in principle nonlinear, function of the components x and u that reflects the modality of interactions between the controller and the existing species. The equation for the dynamical evolution of the controller u displays two distinct contributions. The first represents a self-reaction term, assumed to be linear just for ease of presentation. The nonlinear self-dynamics of the controller u can be readily considered, with no further technical complication. The rate of change of u is assumed to be contextually driven by a global forcing that senses the relative distance of x from its deputed equilibrium . The parameters = (β1, β2, …, β) and ρ will prove central in enforcing the stabilization of the prescribed fixed point. A few comments are mandatory to fully appreciate the generality of the proposed framework, beyond the specific choices made for purely demonstrative purposes. Let us begin by remarking that the controller u can represent an artificially engineered component or, equivalently, belong to an extended pool of interacting populations. In the scheme here imagined, it is assumed that the values of u and x ∀i, are accessible to direct measurements at any time and that this information can be processed to set the controller dynamics. This is largely reasonable for experiments that run under protected conditions like, e.g., the study of microbial dynamics in laboratory reactors, but certainly less realistic for applications that aim at in vivo multidimensional systems, think for instance to genetic regulatory circuits. The dynamical equation for u can, however, be amended to a large extent and with a great deal of flexibility, depending on the target application and the structural specificity of the employed controller, while still allowing for an analogous methodological treatment. As a matter of fact, we can equivalently assume a generalized equation for the controller of the type where f(u*) = 0 and . The dynamics of the original, unsolicited, components and the functional form that specifies the controller feedback bear unequivocal universality traits [15].

The global fixed point (*, u*) of the controlled system (2) should match the following constraints which, provided the and u* are assigned, ultimately set the values of the parameters α. Conversely, as we shall illustrate in the following, one could assume the parameters α as a priori known and infer via Eq (3) the fixed point(s) to be eventually stabilized. The next step in the analysis aims at ensuring the stability of the selected fixed point. This will be achieved by acting on the residual free parameters and ρ. As routinely done, we perturb the equilibrium solution as , u = u* + w and Taylor expand Eq (2) assuming the imposed disturbances small in magnitude. At the linear order of approximation one obtains: where is a N-dimensional column vector of components . The N × N matrix is defined as:

The fixed point (*, u*) is linearly stable if all eigenvalues of the Jacobian matrix have negative real parts. The associated characteristic polynomial can be cast in the equivalent, affine in the ρ-parameter, form: that is reminescent of the celebrated root locus method [32]. Here, [adj()] = (−1) det[()(] denotes the adjugate of matrix , ()( being the minor of obtained by removing the i-th row and the j-th column. The polynomials and have respectively degrees N+1 and N − 1, without loss of generality can be chosen so that assumes this form. With a slight abuse of language we will refer to as poles the roots p of the polynomial and zeros the roots z of . Notice that for ρ = 0 the eigenvalues of the Jacobian correspond to the N+1 poles p. These latter quantities are uniquely determined, once the fixed point (*, u*) has been assigned. In particular it cannot a priori be ensured that the real parts of all p are negative, as stability would require. In other words, when ρ = 0, we can enforce the desired fixed point into the system but cannot guarantee its stability. On the other hand, for ρ → ±∞, N − 1 eigenvalues of tend to the zeros z, which depend self-consistently on the free parameters . As we shall show hereafter, it is in principle possible to assign the β to force the real parts of all z to be negative. The two remaining eigenvalues of matrix , in the limit of large |ρ|, diverge to infinity in the complex plane. More precisely, they travel along opposite directions following a vertical (resp. horizontal) asymptote, if ρ is bound to the positive (resp. negative) semiaxis. To confer stability in the limiting case ρ → ∞ where N − 1 eigenvalues of coincide with the roots of , it is therefore sufficient to (i) operate a supervised choice of and (ii) impose the condition that yields a vertical asymptote (ρ → +∞), while, at the same time, requiring that this latter intersects the negative side of the real axis. In this respect, it is important to remark that the intersection occurs in the point of abscissa . Hence, the idea is to interpolate between the two limiting cases ρ = 0 and ρ → ∞ by determining the minimal value ρ of ρ beyond which the desired fixed point becomes stable. The existence of the threshold ρ that makes the imposed fixed point attractive for any ρ > ρ is obvious, being stability already assured in the limiting setting ρ → +∞. For the sake of clarity we reiterate that this amounts to selecting Re(z) < 0 for all k and further imposing λ0 < 0, by properly assigning the free parameters .

In order to study the assignability of the zeros z by means of , let us recall that for a generic square matrix , where c stands for the coefficients of the characteristic polynomial of , namely . The polynomial can be consequently written as: It is hence straightforward to establish a direct relation between the parameters and the vector of coefficients : that can also be equivalently stated as: where is the matrix defined by: The suited vector is thus obtained from Eq (7), provided matrix is invertible. This latter request defines the condition of controllability for the scheme that we have implemented (see Supporting Information, SI, for a discussion that aims at positioning this observation in the context of standard control theory [33]). Notice that for obvious consistency reasons must have real entries. This follows naturally if one chooses the zeros z to be real or complex conjugate in pairs, which implies that the coefficients d of the polynomial (see Eq (5)) are real. All other quantities involved are real by definition.

Summing up, the devised strategy consists of the following steps. First, the fixed point is selected and the parameters frozen to their respective values as specified by Eq (3). Then the complex roots z are chosen so that Re(z) < 0 for all k while, at the same time, matching the condition that makes the vertical asymptote cross the horizontal axis with a negative intercept. As we will clarify when discussing the applications, the z can be chosen to coincide with the poles p, except for punctual modifications whenever Re(p) > 0. Notice however that z should be real or come in conjugate pairs, as the coefficients d are, by definition, real. Once the roots z have been fixed, one can readily compute the associated polynomial coefficients d, and hence proceed with the determination of via Eq (7), provided that the controllability condition holds. Finally, by selecting ρ > ρ > 0 we obtain a linearly stable fixed point (*, u*) for the controlled dynamics Eq (2).

Results

Testing the control method: From synthetic gene network to real microbiota dataset

As a first application of the above technique, we will study the dynamics of an artificial gene network [34-37]. In our example the network of connections is a regular tree with branching ratio r = 4. It is further assumed that the genetic activation between nodes i and j is described in terms of a Hill function, with cooperation coefficient equal to 2. In formulae, A = 1 and . Negative regulation loops are also accommodated for. These latter could, in principle, be modeled by assuming paired interactions of the type 1 − g(x), while still setting to one the relative entry of the connection matrix. As described in the SI, we can equivalently set A = −1, while assuming interactions to be modulated by g(x) as indicated above. At the same time, the reaction part should be modified with an additional term, η, counting the number of negative loops that affects node i. More specifically, f(x) = −γ x+η, where the first term mimics constitutive degradation. In our tests, matrix contains an identical number of randomly assigned ±1. The parameters γ are random variables uniformly distributed over the interval [0, 1]. Working in this setting our aim is to control the equilibrum state of the system and thus shape the pattern of asymptotic activity. For this initial application we choose to operate with a simple linear control and set h(x, u)≡h(u) = u, for all i. In this case, u could e.g. represent the density of a suitable retroviral vector used to infect specific cell lines [38]. To provide an immediate graphical illustration of the power of the method, we set to stabilize two distinct fixed points. In the first example (see Fig 1, upper left panel) the control is designated so as to enhance the degree of activity of the peripheral nodes of the tree. These latter are characterized by a similar value of the activity, apart for slight randomly superposed fluctuations. Similarly, the nodes that define the bulk of the tree display a shared degree (except for tiny stochastic modulation) of residual activity. In the upper right panel of Fig 1, the dual pattern is instead obtained and stabilized: the peripheral nodes are now being silenced and the activity concerns the nodes that fall in the center of the tree. Lower panel of Fig 1 displays the root locus diagram relative to the situation reported in the second panel. By properly tuning ρ above a critical threshold ρ, we can enforce the stability of the obtained fixed point. Two eigenvalues diverge to ±∞ following a vertical asymptote in the complex plane. For each chosen fixed point that is being stabilized the zeros z can be selected so as to make the asymptote intercept the horizontal axis in the left-half of the plane.

Fig 1

Stabilization of different fixed points for an artificial gene regulatory network.

Upper left panel: the control is modulated so as to enhance the activity of the peripheral nodes of the tree, as compared to the inner ones. Upper right panel: the control makes now the bulk nodes more active as compared to the peripheral ones. Lower panel: the root locus diagram relative to the situation displayed in the second panel is plotted. Blue circles stand for the position of the complex eigenvalues when ρ = 0, while green crosses identify the eigenvalue obtained for ρ → ∞. The vertical red line represents the asymptote that attracts two of the modified eigenvalues, when ρ → ∞. The red dots show the computed spectrum, calculated when increasing ρ. In this case the matrix contains an identical number of ±1 entries. These are randomly assigned and kept unchanged for all tests performed. The figure on the right is a zoom of the plot displayed on the left.

As a second application of the proposed control strategy, we set to study the dynamics of the gut microbiota [6]. The intestinal microbiota is a microbial ecosystem of paramount importance to human health [39]. Efforts are currently aimed at understanding the microbiota ability to resist to enteric pathogens and assess the response to antibiotics cure of intestinal infections. Recent advances in DNA sequencing and metagenomics make it possible to quantitatively characterize the networks of interactions that rule the dynamics of the microbiota ecosystem. This was for instance achieved in [40] by analyzing available data on mice [41] with an innovative approach which combines classical Lotka-Volterra model and regression techniques. Eleven species were identified and thoroughly analyzed in terms of self and mutual dynamics.

In the following we shall apply the method here developed to control the dynamics of the whole microbioma [40] or a limited sub-portion of it. In this specific case, the self-dynamics is assumed to be logistic, namely f(x) = x(r − s x), while g(x, x) = x x. The constants r and s are provided in [40] and follow from direct measurements. The weighted matrix of connections presents both positive and negative entries, assigned according to [40]. Finally, h(x, u) = ux. The results of the analysis are organized under different headings that reflect the three distinct control strategies explored.

Stabilizing an unstable fixed point by means of an external controller (Case A)

Consider the system of 11 species, as defined in [40] (see SI for a discussion on the bacterial species involved). For illustrative purposes, we will restrict the analysis to all sub-systems that combine 5 out of the 11 species analyzed in [40]. The fixed points for the obtained 5 species systems are calculated. Those displaying positive concentrations are then retained for subsequent analysis. The stability of each selected fixed point is established upon evaluation of the spectrum of the Jacobian of the reduced dynamics. In Fig 2(a) the histogram of (λ), the largest real parts of the recorded eigenvalues, is plotted: several fixed points exist that correspond to unstable equilibria. Starting from this setting, we will introduce a suitably shaped controller, following the above discussed guidelines, in order to stabilize a slightly perturbed version of an originally unstable fixed point, see pie charts in Fig 2(a). Denote by the fixed point to be eventually stabilized and consequently assign the parameters α so as to match Eq (3). The spectrum of the Jacobian matrix obtained for ρ = 0 (blue circles in Fig 2) protrudes into the right half-plane. More specifically, one eigenvalue exhibits a positive real part, so flagging the instability that one aims to control. At variance, the crosses in Fig 2(b) stand for the roots z of , and fall in the left side of the complex plane. The vertical (red, in Fig 2(b)) line identifies the location of the two residual eigenvalues of the Jacobian matrix, when ρ → ∞. By tuning the parameter ρ, one can continuously bridge the two above limiting settings, as graphically illustrated in Fig 2(b). When ρ > ρ ≃ 0.01, the eigenvalues populate the left half-hand plane and stability is, therefore, gained.

Fig 2

Control strategies applied to the microbiota dataset.

Panel (a): a reduced 5 species subsystem of the microbiota is considered (case A) and all possible fixed points computed. Only those displaying non-negative concentrations are retained and their stability assessed. In the main figure, the histogram of (λ), the largest real parts of the eigenvalues obtained after the linear stability analysis, is depicted. The two pie charts refer to the initially unstable fixed point (upper chart) and the stabilized equilibrium (lower chart). Panel (b): the root loci diagram relative to the case discussed in panel (a), is shown. Blue circles identify the position of the complex eigenvalues when ρ = 0, while green crosses stand for the eigenvalues obtained in the limit ρ → ∞. The vertical red line is the asymptote that eventually attracts the two residual eigenvalues. The red dots show the computed spectrum, when progressively increasing ρ. Panel (c): the goal is here to reduce the concentration of the pathogen species, C. difficile, by employing as controller one of the species that compose the microbioma (case B). The concentration of C. difficile is monitored over time for three different control strategies, turning on the control at the same time (t = 20 days). The insertion of the species of uncl. Lachnospiraceae provoques a substantial reduction (50%) of the pathogen concentration, as also displayed by the enclosed pie charts (for interpreting the color-code refer to panel (a)). Panel (d): we now modify a stable fixed point, by driving to extinction one of the existing population, the pathogen C. difficile (here species 6), with an indirect control strategy (case C). The obtained concentrations are reported in the left graph (pluses) and confronted with the initial unperturbed solution (diamonds). As anticipated . The components of are plotted in the right graph. Notice in particular that α6 = 0.

Acting with one species of the pool to damp the concentration of the pathogens (Case B)

Select now a stable fixed point, mixture of five distinct species. One of them is Clostridium difficile, a species of Gram-positive spore-forming bacteria that may opportunistically dominate the gut flora, as an adverse effect of antibiotic therapy. As controller we shall here employ one of the other 6 species that compose the microbioma [42, 43]. The aim is to drive the system towards another equilibrium, stable to linear perturbations, which displays a decreased pathogen concentration. In this case the parameters are determined a priori, once the control species has been identified. Denote by the reduced 5 × 5 matrix that specifies all paired interactions between the pool of populations involved in the initial fixed point. The equilibrium solution that can be attained by the controlled system is determined as , and clearly depends on the species used as controller. The only meaningful solutions are those displaying non negative components . In the example depicted in Fig 2(c) only three solutions can be retained, namely the ones obtained by using uncl. Lachnospiraceae, uncl. Mollicutes and Enterococcus as respective control. In one of the inspected cases, the amount of C. difficile is found to reduce, when the control is turned on. The asymptotic concentration that is eventually attained is sensibly lower than the one initially displayed. The pie charts in Fig 2(c) represent, respectively, the initial fixed point and the final stationary equilibrium, as shaped by the control in the most beneficial case, i.e., when the concentration of C. difficile is seen to shrink. The root locus plot obtained for this specific case is reported in the SI. Importantly, the discussed scheme can be straightforwardly modified so as to account for a generic nonlinear self-reaction dynamics for the control species, e.g., a logistic growth, that could replace the linear Hookean-like term assumed in Eq (2).

Driving to extinction one species, the other being the target of the control (Case C)

As an additional example, we wish to modify a stable fixed point of the dynamics, by silencing one of the existing populations with an indirect control. In other words we shall introduce and stabilize a novel fixed point, that displays a negligible residual concentration of the undesired species, by acting on the other species of the collection. This is for instance relevant when aiming at, e.g., eradicating a harmful infection that proves resistant to direct therapy. With this in mind, we consider a reduced ecosystem consisting of 6 species, selected among the 11 that define the microbiota. A stable fixed point exists (black diamonds in Fig 2(d), left panel) which displays a significant concentration of C. difficile, the pathogen species. Assign to this latter species the index 6. We now insert a controller which cannot directly interfere with C. difficile. This amounts, in turn, to setting to zero the corresponding component of vector (α6 = 0). We then require the concentration of the C. difficile to be small, i.e., . This latter condition translates into a constraint that should be matched by the other 5 species, namely . Given , the components α, with k ≠ 6, are chosen so as to match the constraint . A possible solution of the problem is reported in Fig 2(d): in the left panel (plus symbols) the components of the fixed point stabilized by the control are shown. As anticipated, the concentration of C. difficile is small. The right panel of Fig 2(d) shows the components of the vector that specify the characteristics of the introduced controller. Notice that α6 = 0 so that the controller is not directly influencing the rate of production of C. difficile.

Discussion

We would like to draw the attention on the interpretation of . As stated earlier, characterizes the strength of the coupling between the controller u and every single species of the system to be controlled. An alternative interpretation is however possible: u could represent a mixture of different species and the components of incorporate the relative abundance of the mixed compounds. In light of the above, also the previously discussed control schemes which apparently assumed dealing with an artificially designed control, could be realized via a proper mixture of exisiting microbiota species so as to achieve the coupling corresponding to the desired fixed point.

Notice also that the control scheme here developed could be in principle exploited to drive the system towards a stable fixed point of the unperturbed dynamics, starting from out-of-equilibrium initial conditions. To achieve this goal u* needs to be set to zero, thus requiring that the controller is turned off at equilibrium. In this case, and are not subjected to specific constraints, as the existence and stability of the desired equilibrium are a priori granted. Such parameters could hence be chosen so as to reflect the specificity of the target system. In the annexed SI we demonstrate this intriguing possibility.

Conclusion

Summing up, we have here proposed and tested a method to control the dynamics of multidimensional systems on a complex graph. The original system is made up of N interacting populations obeying a set of general equations, which bear attributes of universality. One additional species, here referred to as the controller, is inserted and made interact with the existing constellation of species. By tuning the strength of the couplings (or equivalently the composition of the inserted controller), we can drive the system towards a desired equilibrium. The stability of the achieved solution is enforced by adjusting the parameters that ultimately govern the rate of change of the controller. Methodologically, we make use of the root locus method which can be naturally invoked once the control problem is suitably formulated. The tests that we have performed, both synthetic and drawn from real life applications, demonstrate the versatility and robustness of the proposed scheme. This latter configures therefore as a viable and innovative tool to tackle a large plethora of inter-disciplinary systems, from life science to man-made applications, that should be stably driven towards a desired configuration. In this current implementation, and for purely pedagogical reasons, the control assumes that the state of the system is accessible to direct measurement. Relaxing this working hypothesis is a possibility that we shall explore in a future contribution.

Control of multidimensional systems on complex network: Supplementary information.

(PDF)

Click here for additional data file.