Spectral control for ecological stability
Giulia Cencetti Franco Bagnoli Giorgio Battistelli Luigi Chisci Duccio Fanelli
cond-mat.stat-mech physics.soc-ph q-bio.PE

	
	
	
	
	 
	Dipartimento di Ingegneria dell'Informazione, Università di Firenze,
		Via S. Marta 3, 50139 Florence, Italy
	INFN Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italia
	  Università degli Studi di Firenze, Dipartimento di Fisica e Astronomia and CSDC, via G. Sansone 1, 50019 Sesto Fiorentino, Italia
	INFN Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italia
	 
	Dipartimento di Ingegneria dell'Informazione, Università di Firenze,
		Via S. Marta 3, 50139 Florence, Italy
	 
	Dipartimento di Ingegneria dell'Informazione, Università di Firenze,
		Via S. Marta 3, 50139 Florence, Italy
	 Università degli Studi di Firenze, Dipartimento di Fisica e Astronomia and CSDC, via G. Sansone 1, 50019 Sesto Fiorentino, Italia
	INFN Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italia
	
	
	
	
			A system made up of N interacting species is considered.  Self-reaction terms are assumed of the logistic type. Pairwise interactions take place among species according to different modalities, thus yielding a complex asymmetric disordered graph. A mathematical procedure is introduced and tested to stabilise the ecosystem via an ad hoc rewiring of the underlying couplings. The method implements minimal modifications to the spectrum of the Jacobian matrix which sets the stability of the fixed point and traces these changes back to species-species interactions. Resilience of the equilibrium state appear to be favoured by predator-prey interactions.
			
			
	

	Spectral control for ecological stability
    Duccio Fanelli
    December 30, 2023
=========================================





§ INTRODUCTION





A vast plethora of physical phenomena occurring in nature are due to collective dynamics, which spontaneously emerge at the macroscopic level in systems made up of microscopically interacting constituents. When the entities that compose the whole set are subject to specific self-reactions and, at the same time, diffuse across the embedding spatial medium, the system goes under the name of reaction-diffusion  <cit.>. Examples are invasion models in ecology <cit.>, epidemic spreading <cit.>, and also the celebrated Turing patterns that arise, for instance, from the dynamical interplay between reaction and diffusion in a chemical system <cit.>.


Usually, reaction-diffusion models are defined on a regular lattice, either continuous or discrete <cit.>. In many cases of interest, it is however necessary to consider the units moving on a complex network <cit.>, where each node represents a physical location; individual elements are constantly displaced from one node to another following diffusion rules. Spatial ecological systems can be schematized by resorting to a network that exemplifies existing routes linking different habitats, infection spreading customarily requires accounting for a transportation network, and also chemical reactors and biological cells can be considered as coupled through complex networks.



A reaction-diffusion model defined on a complex network can display a stationary stable equilibrium. Stability is of paramount importance as it relates to resilience, the ability of the system to counteract external perturbations that would tend to get away from the existing equilibrium. It is therefore crucial to devise possible strategies aimed at interfering with the system of interest so as to enforce the desired stability. This can be achieved by acting on the local dynamics or, more interestingly for what it follows, by reshaping the underlying network of spatial connections. This latter possibility has been addressed in <cit.>, via a minimally invasive procedure targeted to modifying the spectrum of the Laplacian operator which governs the diffusive, hence linear, couplings. This, in turn, implies rewiring the links, and consequently recalibrating the weights, of the underlying graph of connections.
 

Dynamics on networks is however central to other realms of investigation, also when the system being examined is not made spatially extended. In this case, long-ranged (non local) interactions, as encoded in the structure of the assigned network, might follow non-diffusive rules.  In ecosystems inter-species interactions are assumed to be 
mediated by pairwise, hence quadratic, exchange, to some level of approximation. Each population (species) is characterised by a self-reaction dynamics, typically described via a suitable nonlinear function of the concentration amounts, and different populations can be abstractly assigned to given nodes of a virtual graph. The network then represents the interactions between different species and the sign of the weighted entries of the associated adjacency matrix define the nature of coupling (competitive, cooperative, predator-prey, etc...).

The stability of the ecosystem, or its resilience, is essential in many respects. In the past decades many efforts have been made to understand the principles that rule the stability of a complex ecosystem. These concepts have been originally addressed by Robert May <cit.> in a seminal paper that paved the way to a completely new field of investigation, still very productive and fertile. These studies resulted in multiple attempts of providing methods to make an ecological community structurally stable <cit.>. May <cit.> analysed a system described by N variables (N interacting species) obeying a set of differential equations. The stability analysis is performed by linearising the equations in the neighbourhood of an equilibrium point, whose stability depends on the spectrum of the interaction matrix. May's analysis focuses on this latter, eventually bearing to the challenging statement that, in short, the more complex the more unstable is the system. Recent work by Allesina et al. <cit.> provides an implementation of May's ecosystem accounting for well defined (non random) interactions. These are competitive, mutualistic and predator-prey and according to Allesina et al. play different roles in the stability of the ecosystem. Remarkably, the presence of predator-prey relations has a stabilizing effect. In Coyte et al. <cit.> stability of a microbioma ecosystem is obtained by allowing for sufficiently weak couplings. 


Starting from these premises, in this paper we provide an alternative approach to ecological stability by developing a self-consistent mathematical strategy which implements a spectral control algorithm. The method builds on the technique developed in <cit.> and extends its domain of applicability, beyond diffusion mediated (linear) processes to the interesting setting where pairwise, hence non-linear, non-local interactions are considered. In doing so, we will contribute to identifying the key topological features that should be possessed by a stable (resilient) ecological network. To anticipate our findings we will show that predator-prey interactions exert a beneficial role in terms of stability in qualitative agreement with the results reported in  <cit.>. Further weak interactions tend to favour the overall stabilization, as observed in <cit.>.

The paper is organized as follows. In Section <ref> we will introduce the model and define the reference mathematical setting. In Section <ref> the control scheme is discussed in detail and we will also show how the spectral modifications trace back to actual changes in the matrix of interactions. In Section <ref> the method is applied to a simple bidimensional Lotka-Volterra model, while the extension to arbitrarily large systems is discussed in Section <ref>. Finally, in Section <ref> we sum up and draw our conclusions. Relevant technical material is provided in Appendices. 






§ THE MODEL



        We shall hereafter consider the coupled evolution of N species and denote by x_i, i=1, ..., N, their associated concentrations. We will operate under the deterministic viewpoint and deliberately omit any source of stochastic disturbance, be it endogenous (demographic noise) or exogenous (external perturbation). The evolution of the ecosystem is hence described by the following set of first-order differential equations: 
		
	
    ẋ_i=x_i(r_i-s_ix_i) + x_i∑_j≠ iA_ijx_j   i,j=1,...,N.

	
The self-reaction term is here assumed logistic, for pedagogical reasons. We will subsequently relax this working hypothesis and generalize the analysis so as to account for an extended family of nonlinear reaction terms. In the above equations, r_i stands for the intrinsic growth rate of species i, while s_i is inversely proportional to the assigned carrying capacity. Matrix A, in general asymmetric, defines the interactions among species. More specifically, the scalar entry A_ij encodes the effect exerted by species j on species i. The magnitude of A_ij weights the strength of the interactions. The sign of A_ij, respectively A_ji, defines the specific nature of the interaction
 between species i and j. Adopting a wording which is inspired to ecological applications: exploitation (+,-), competition (-,-), cooperation (+,+), commensalism (+,0), amensalism (-,0) or null interaction (0,0). The coupling among species is shaped by a quadratic term, which scales like the product of relative concentrations.  This implies assuming the interaction to be mediated by pair exchanges, as it is customarily the case in ecology.  On a more fundamental level, the coupling term here introduced will enable us to generalize the analysis reported in <cit.> beyond standard diffusion. 
 
 
  Let us denote by x^* the fixed point of the dynamics and assume for the sake of simplicity that all entries x^*_i are different from zero. In formulae we have:	
	
	
    r_i-s_ix_i^*+∑_j≠ iA_ijx_j^*=0   i=1,...,N
	
which essentially defines the sole non trivial equilibrium eventually attained by the system under scrutiny. As alluded to in the introductory Section, we are here concerned with the resilience of the system, namely its inherent capability to recover from perturbations. Stated it differently, we shall elaborate on the conditions which make the 
fixed point stable, according to a linear stability analysis. To this end, it is convenient to introduce the rescaled variable y_i≡ x_i/x_i^*. In the new variables the fixed point is homogeneous and reads 
y_i^*=1 ∀ i. The dynamics of system (<ref>) can be cast in the form:
 
	
    ẏ_i=y_i(r_i-s̃_iy_i+∑_j≠ iB_ijy_j)

where s̃_i=s_ix_i^*, B_ij=A_ijx_j^* and the equation (<ref>) assumes the form r_i-s̃_i+∑_j≠ iB_ij=0.

To assess the stability of the fixed point, we set y_i=1+v_i and Taylor expand at the first order in the perturbation amount. This yields the following linear system for the evolution of the imposed perturbation:

	
    v̇_i=-s̃_iv_i+∑_j≠ iB_ijv_j ≡∑_jC _ijv_j.


The stability of the fixed point is ultimately controlled by the eigenvalues of the matrix C obtained by adding the elements -s̃_i on the diagonal of matrix B. The system is unstable when at least one eigenvalue of C has a positive real part. In the following, we will discuss a control procedure to stabilise a fixed point, that is initially engineered to be unstable. The method builds on the technique discussed in <cit.> and aims at reshaping the coupling among interacting species. As a side observation, which will become crucial in the forthcoming analysis, we notice that the fixed point condition (y_i^*=1 ∀ i) translates into a constraint for C: summing the elements of C relative to row i one should recover -r_i, i.e. ∑_j C_ij = -r_i.


		
	
	

§ TOPOLOGICAL CONTROL SCHEME

	
	
	
The proposed control strategy aims at modifying an initially unstable system of the type described above to yield an equivalent analogue which preserves the  form (<ref>) while admitting a stable non trivial equilibrium. As it shall be argued hereafter, we can either enforce the stability of the original, assumed unstable, fixed point, or, alternatively, steer the system towards a different equilibrium. In the former case we shall also alter the original carrying capacity (a parameter which indirectly encodes for the interaction with the surrounding environment), while in the latter the envisaged protocol will solely impact the network of interspecies couplings, leaving unchanged individual reaction parameters. In both cases,  the intrinsic growth rates r_i, are kept unvaried. Stable ecological networks display remarkable topological characteristics, as discussed in a series of  papers devoted to this topic <cit.>. Operating along this line, we will isolate and discuss a selected gallery of features that appear to be recurrently shared by the ecological networks stabilized as outlined in the following. 

Since, by definition, the structure of (<ref>) is invariant under the foreseen procedure,  inspecting the linear stability of the ensuing equilibrium implies dealing with a system  of the type v̇_i=∑_jC'_ijv_j, where C' is obtained from C, defined as in equation (<ref>), via the devised control algorithm. Requiring the sought stability is, in turn, equivalent to constrain the spectrum of C' in the left portion of the complex plane, namely to set the real part of the associated eigenvalues to negative values. Our goal, pursued hereafter, is to elaborate on a rigorous mathematical procedure, which is both anchored to first principles and potentially minimally invasive, to derive C' from C. Importantly, the obtained matrix C' should match the condition ∑_j C'_ij = -r_i, for the homogeneous fixed point to exist in terms of the rescaled variables y_i (recall that, by hypothesis, r_i is  frozen to its original value). The effect of the control will be then gauged by tracing the modifications back to the underlying nonlinear framework, i.e. by evaluating the impact produced on the relevant dynamical parameters. 

As a  first step in the analysis, we write the linear equation (<ref>) in the equivalent form:

      
    v̇_i=-r_i v_i+∑_jD_ijv_j.

where the definition of D follows trivially. Recalling that ∑_j C_ij = -r_i by virtue of the aforementioned fixed point condition, it is immediate to conclude that 
D is a zero-row-sum matrix, namely ∑_j D_ij =0. The next step is to diagonalize matrix D. Formally, we set 
Φ^-1 DΦ=Λ, where Φ is the matrix whose columns are the eigenvectors of D and Λ the diagonal matrix formed by the corresponding eigenvalues.  Diagonalizability of matrix D is hence a necessary requirement for the method to hold.
The idea is now to calculate the (minimal) shifts δΛ^(α), α=1,...,N, to be applied to the eigenvalues of  Λ^(α) of matrix D, for the homogeneous fixed point  y_i^*=1 ∀ i to prove linearly stable. Recall that the stability of this latter fixed point is eventually dictated by the spectrum of C (or, more precisely, by its controlled version C'), from which the zero-row-sum counterpart D originates. The needed corrections δΛ^(α) are organized in a N × N diagonal matrix δΛ (with δΛ_kk=δΛ^(k)) to be added to matrix Λ. The key point is how to choose the 
entries of  δΛ for the control to return an effective, moderately intrusive (in terms of the modifications apported on the spectrum),  stabilization. 

To answer this question, we proceed as if the original eigenvalues Λ were perturbed by the finite amount δΛ and recover a matrix D', which displays the modified spectrum, via the inverse transformation D'≡Φ(Λ+δΛ)Φ^-1. By construction D' commutes with D, as the two matrices share the same set of eigenvectors. Notably, the corrections δΛ can be chosen such that matrix D' is also zero-row-sum, as D is. This is rigorously proven in Appendix <ref>, building on the derivation reported in <cit.>. Moreover,  D' is real: this property is also inherited from D, as shown again in Appendix <ref>. We are thus brought back to the linear problem:

      
    v̇_i=-r_iv_i+∑_j D'_ijv_j

or, equivalently to v̇_i=∑_jC'_ijv_j, where:
	
    C'_ij= {[ D'_ij      if i≠ j; D'_ii-r_i  if i=j.;                    ].


By construction ∑_j C'_ij = -r_i, since  ∑_j D'_ij = ∑_j D_ij  = 0. This is a crucial observation which makes it possible to interpret C' as the Jacobian matrix 
associated to a rescaled nonlinear problem of the type (<ref>) where parameters r remain unchanged. 	The zero-row-sum-property of matrix D' is not a necessary condition to obtain stability but represents a useful requirement to help interpreting the results in terms of the original variables.
 
We now return to discussing the selection of the elements of the shift matrix δΛ. These latter are to be chosen so as to constrain the spectrum of C' to the left hand side of the imaginary plane, thus ensuing the desired stability.  
From (<ref>), it is clear that the eigenvalues of D' are positioned, in the complex plane, on the right of  those stemming from matrix C'. The relative separation between the two respective spectra can be somehow quantified through r. To make this observation rigorous, we recall the celebrated Gershgorin theorem <cit.>:
the eigenvalues of a given matrix are included in disks defined by the elements of the matrix itself. More specifically, the i-th Gershgorin disk of matrix D' corresponds to the relative disk of matrix C', translated to the right by the scalar quantity r_i. Unfortunately, it is not trivial to relate the index i (running on the nodes) to the eigenvalues (sorted with the index α).  To enforce stability, and assuming the worst case scenario, we shall assign the (real) shifts δΛ^(α) so that
all eigenvalues of D' have their real part smaller than r_min, the minimum of all the entries of vector r. In practical terms, the imposed corrections δΛ^(α) are chosen as:

	
    δΛ^(α)= {[ R - Λ^(α)_Re   if Λ^(α)_Re>r_min;                       0otherwise;                                  ].

where the scalar quantity R has been introduced such that R<r_min. To help visualizing the whole procedure we report an illustrative example in Figure <ref>, without insisting on the specific selection of the involved parameters. The eigenvalues of the original matrix C are displayed in the complex plane with (blue) plus symbols: the homogeneous fixed point of the rescaled equations (<ref>) is therefore unstable, as the spectrum protrudes in the right half-plane. The eigenvalues of  D are shown with (yellow) triangles and extends on the right of the vertical dashed line, which is traced at r_min. The (purple) circles stands for the eigenvalues of the controlled matrix D': as anticipated, they are confined on the left of the vertical dashed line, the closer to the line the less invasive the control imposed on the population of unstable modes.  Finally, and as predicted, the spectrum of  the controlled Jacobian C' is contained in the negative half-plane, thus implying asymptotic stability. 




As it should be clear from the above,  the control protocol assumes dealing with constant  r_i parameters. Given this constraint, two viable strategies are envisaged to re-parametrize the original system, in light of the outcome of the control scheme. We have in fact 

    		
    D'_ij= {[       r_i-s_i'x_i^*'if i=j; B'_ij = A'_ijx_j^*'if i≠ j;                            ].
 
which allows in principle to define the new coupling strengths, as encoded in A', the novel fixed point x^*' and the modified inverse carrying capacities s'. A first strategy to finalize the transformation suggests leaving the parameters s unchanged, namely s'= s. In practical terms, we assume that the reaction parameters, which characterize the dynamics of each species when evolved on a isolated patch, remain unchanged. The ecosystem can achieve stabilization, by just reshaping the underlying networks of interlaced dependencies. From equation (<ref>), we have therefore:

	
    [                   x_i^*'=r_i-D'_ii/s_i; A'_ij=D'_ij/x_j^*'=D'_ijs_j/r_j-D'_jj.;                                        ]

	
Setting paired interactions as specified by matrix A'  guarantees the stability of the associated, and consistently modified, fixed point x^*'.  The quantities x_i^*' should be positive defined (at least when ecological applications are concerned), which in turn translates into the additional requirement  

    r_i-D'_ii>0 ∀ i.

	
The second strategy consists of modifying the parameters s, together with the matrix A, leaving unchanged the fixed point x^*. The carrying capacity s is prone to environmental influences, and, as such, it can be imagined to be tunable with some degree of realism. This is opposed to the growth parameter r,  constrained, among other factors, by  species genetics, and thus assumed constant throughout the procedure. From equation (<ref>):
	
    [ s_i'=r_i-D'_ii/x^*_i;   A'_ij=D'_ij/x_j^*. ]
 
	The additional condition  s_i'>0 should be imposed, which again amounts to requiring  equation (<ref>) to hold. In other words, condition (<ref>) is  a general    constraint that the control scheme is bound to verify, for the specific ecologically inspired application, here discussed. We shall refer to condition (<ref>) as to the applicability constraint and elaborate on its implications hereafter. The idea is to find a suitable value for R in order to match the condition (<ref>), and, hence, to make the control scheme applicable. To begin, let us assume that we are allowed to modify all eigenvalues of the considered spectrum, and not just the limited sub-set that triggers the system unstable.  Then, it is enough to impose δΛ^(i)≤s̃_min ∀ i, which is in principle always possible, in order to automatically verify (<ref>). The system is therefore always controllable when all eigenvalues are to be modified.


Consider now the more interesting case where a subset ℳ of elements of the whole spectrum is the target of the control. One has to face the following restrictions:
    
        
  * if an index i exists such that ∑_j∈ℳΦ_ijΦ_ji^-1<0 and s̃_i+ ∑_j∈ℳΦ_ijΛ_Re^(j)Φ_ji^-1 - r_min∑_j∈ℳΦ_ijΦ_ji^-1<0, then constraint (<ref>) is never matched and the system is not controllable with the  above discussed technique.
        
        
  * if an index i exists such that ∑_j∈ℳΦ_ijΦ_ji^-1>0 and s̃_i+ ∑_j∈ℳΦ_ijΛ_Re^(j)Φ_ji^-1 - r_min∑_j∈ℳΦ_ijΦ_ji^-1<0 (suppose node indices are sorted so that such an index i is in the subset of nodes 1,...,ñ with ñ < N), than the applicability condition (<ref>) is verified only if the following statement holds true:
        
    max_i∈[1,ñ](r_min - s̃_i+ ∑_j∈ℳΦ_ijΛ_Re^(j)Φ_ji^-1/∑_j∈ℳΦ_ijΦ_ji^-1)
                ≤min_i∈[n+1,N](r_min - s̃_i+ ∑_j∈ℳΦ_ijΛ_Re^(j)Φ_ji^-1/∑_j∈ℳΦ_ijΦ_ji^-1)

        where [n+1,N] numbers the set of indices i for which ∑_j∈ℳΦ_ijΦ_ji^-1<0.
    

The above conditions (derived in Appendix <ref>) are to be carefully checked before attempting to control the system under scrutiny via the procedure that we have here illustrated and which is ultimately aimed at recalibrating the weights of the underlying couplings.
	
	
	
Before concluding this Section, we briefly mention an alternative control strategy which builds on the already mentioned  Gershgorin theorem. In the above analysis we have reshaped the networks of contacts and altered either the fixed point x^* or the carrying capacity s_i, while preserving the values of the growth factors r_i. The alternative route to stabilization that we shall hereafter discuss follows a dual path: the only parameters to be tuned are the growth rates r_i. Recalling that C' is real (from definition of (<ref>), being D' real), it is clear that each eigenvalue has its real part smaller than C'_ii+∑_j≠ i|C'_ij|, which is the rigthmost point in the complex plane of the i-th Gershgorin disk for matrix C'. To enforce stability, we could then require that all the Gershgorin circles are included in the left half-plane:
	
    C'_ii+∑_j≠ i|C'_ij|≤0 ⇒ D'_ii-r_i+∑_j≠ i|D'_ij|≤0 ⇒ D'_ii+∑_j≠ i|D'_ij|≤ r_i

	which implies that the Gershgorin disks computed for matrix D' should be contained in the semi-plane constrained, from the right, by the  vertical line  located at  r_i for any selected i, or, better, at the left of r_min, for all i. The shift of the Gershgorin disks can be performed by modifying the diagonal of matrix C, which is equivalent to changing the vector of parameters r. This alternative control strategy proves however more invasive in terms of the perturbation that is produced on the original spectrum. This is clearly testified in Figure <ref>, where the original (unstable by construction) spectrum is compared to those obtained applying the two control strategies outlined above.  

	
	

	
	
In the next Section we will begin by applying the developed method to a simple example, where just two species are made to evolve. This application bears pedagogical interest and it will pave the way to inspecting the general setting, on which we shall report in the subsequent Section. 	

	
	

	

§ STABILIZING THE LOTKA-VOLTERRA DYNAMICS

	
	
With the goal of gaining further insight on the control scheme developed above, we will here consider a simple setting where just two species are made to mutually interact. Hence, we will consider hereafter  N=2 and label with x_1 and x_2 the mean field concentrations of the interacting species. Equations (<ref>) reduce therefore to:

	
    {[ ẋ_1=x_1(r_1-s_1x_1+A_12x_2); ẋ_2=x_2(r_2-s_2x_2+A_21x_1);                             ].

a framework known in the literature as the Lotka-Volterra competition model. To proceed in the analysis we define u_i≡ (s_i / r_i ) x_i, τ≡ r_1t,  and eventually get the following governing equation in the rescaled variables u_1 and u_2:	
		
		
    {[   du_1/dτ=u_1(1-u_1+au_2)≡ f_1(u_1,u_2);                                        ; du_2/dτ=ρ u_2(1-u_2+bu_1)≡ f_2(u_1,u_2);                                         ].

where ρ=r_2/r_1, a≡-A_12r_2/r_1s_2 and b≡-A_21r_1/r_2s_1. It is straightforward to prove that this system admits four equilibria: (u_1, u_2)=(0,0), (u_1, u_2)=(0,1), (u_1, u_2)=(1,0) and (u_1, u_2)=(1-a/1-ab,1-b/1-ab)≡(u_1^*, u_2^*). We shall hereafter refer to the latter equilibrium, the only one to guarantee non trivial asymptotic concentrations for both species. Notice that this is admissible only if u_1^* and u_2^* are positive and finite, which, in turn, implies that the parameter space (a,b) is restricted to the colored region of Figure <ref>, i.e.
	
	
		
  * a>1, b>1
		
  * a<1, b<1, b<1/a
	


	
	
	
	The stability of the selected fixed point is determined by the  Jacobian matrix
	
    J_(u_1^*,u_2^*)=
    	1/1-ab (
    	
    	a-1     a(a-1)
    ρ b(b-1)    ρ(b-1)
     )

	whose eigenvalues are 
	
    λ_±=1/2(1-ab)[(b-1)(1+ρ)±√((a-1)^2(1+ρ)^2-4ρ(1-ab)(a-1)(b-1))].

	The sign of λ_± implies that the fixed point (u_1^*,u_2^*) is unstable for a>1 and b>1 and stable in the complementary domain, as depicted in Figure <ref>.

	
	
	
Working in this simplified setting, it is therefore straightforward to implement, and graphically illustrate, the stabilization protocol, as addressed in the preceding Section. 
Starting from an unstable system corresponds to setting the parameters a and b in the (orange) shaded sub-portion of the first quadrant of Figure <ref>. When implementing the control, panel (b) of Figure <ref>,  the rescaled parameters (a,b) are consequently moved to the other (cyan) shaded domain displayed in Figure <ref> (a), i.e. the region deputed to stability. Different symbols, as depicted in the Figure <ref> (a), refers to distinct choices of the initial model parameters. In all cases, the stabilization is successfully produced and, more importantly, the modified parameters (a,b) tend to cluster in a limited portion of the stability domain, close to the boundary of instability. This observation  suggests again that the devised control acts  by producing a  somehow minimal perturbation to the original model.  Following the alternative control recipe based on the Gershgorin theorem 
produces much more invasive changes (here not shown). 
	
	
The controlled matrix of interaction, as well as the novel set of dynamical parameters, can be readily obtained from the modified quantities a' and b',  following one of the interpretative scenarios discussed with reference to the N-dimensional case.  Both strategies require altering the matrix of couplings A to eventually obtain its modified counterpart here denoted with A'. 
More specifically:


	
		
  * Changing the fixed point x^*.  The 
		 parameters r and s stay unchanged. Hence:
		
		
		
    A'_12   =    -a' r_1 s_2/r_2
    
    		A'_21   =   -b' r_2 s_1/r_1

		

		The new fixed point (x_1^*)' , (x_2^*)' is obtained by solving the self-consistent equations:
		
    r_1-s_1(x_1^*)'+A'_12(x_2^*)' =0
    
    		r_2-s_2(x_2^*)'+A'_21(x_1^*)' =0.

		
		
  * Modifying the parameters s. The fixed point and r are not varied. One gets:
		
		
    {[ a'=-A'_12r_2/r_1s_2'; b'=-A'_21r_1/r_2s_1';                      ].

		which, together with the fixed point condition:
		
    {[  r_1-s_1'x_1^*+A'_12x_2^*=0; r_2-s_2'x_2^*+A'_21x_1^*=0.;                             ].

		allows one to obtain the entries of the matrix A' and the controlled parameters s.
	
	
In the first case we have to make sure that the new fixed point is admissible, which amounts to meeting the conditions (x_1^*)'>0 and (x_2^*)'>0. In the second case, one has to impose s_1'>0 and s_2'>0. It is straightforward to prove that these constraints are always satisfied when a' and b' fall, as they do by definition, in the stability region (1-a'>0, 1-b'>0 and 1-a'b'>0). In the next Section we move on to considering the general setting by working with an arbitrarily large ecological system consisting of N interacting populations.

	
	

§ THE GENERAL CASE

	
	
	Let us now consider the general setting where an arbitrarily large number of  species is made to interact. The initial coupling strengths which exemplify the interaction among distinct populations, as encoded in the elements of matrix A, are initially assigned in such a way that a non trivial fixed point exists and is linearly unstable. This is achieved by generating the random entries A_ij, according to a specific distribution that we discussed in the Appendix and exploiting again the Gershgorin theorem. In the following, we will set N=100 and illustrate the results which are obtained when allowing for the fixed point to be modified by the control procedure. The alternative control strategy, which leaves the original fixed point unchanged, is also analyzed and the results reported in the Appendix. As an interesting outcome, we will show that predator-prey interactions exert a stabilizing effect, as already pointed out in <cit.>. 
	
	In Figure <ref>, the original and modified spectrum are displayed in the complex plane, for one representative realization of matrix A. The symbols are chosen following the same convention adopted in Figure <ref>. The stabilization produced by the control is clearly demonstrated. 
	
	
	
	
	
	To help visualizing how the control shapes the examined ecosystem, we extract from the matrix A' the number of pairs that belong to the different classes, categorized in five six large classes, here recalled for the sake of completeness: exploitation (+,-), competition (-,-), cooperation (+,+), commensalism (+,0), amensalism (-,0) or null interaction (0,0). 
	The frequency of appearance of different classes  is investigated in Figure <ref>. The red plus symbols point to the uncontrolled setting, and reflect the specific rule chosen for generating matrix A and its associated, unstable fixed point (see Appendix <ref>). For the controlled matrix, only three bins are populated, see Figure <ref> (a): interaction modalities that envisage a one directional coupling, or stated differently, a zero entry in matrix A, are absent in the controlled scenario. The topological control activates in fact all pairwise connections, albeit often by a tiny amount.  It is then interesting to silence, a posteriori, in the controlled adjacency matrix A',  the links that are associated to  a weight (in absolute value) smaller than a given cutoff. For a sufficiently small cutoff the stability of the controlled system is preserved: the number of newly added links
	can be hence considerably reduced, by eradicating from the collection those that bear no relevance in light of the modest exerted coupling. The effect of the cut off is visible in Figure <ref> (b): the final distribution of pairs resembles very closely the one generated at the beginning. Remarkably, the number of predator-prey interactions grows at the 
	detriment of the last column of the histogram, implying that the new interactions that are to be established for stability to hold belong to this class, in qualitative agreement with the analysis by Allesina et al. <cit.>.
	
	
		
	
	
	A different view on the effect of the control method can be gained by looking at Figure <ref>. Here, the distribution of the weights associated to predator-prey interactions is plotted before (blue bars and  continuum profile) and after (orange histogram) application of the control. While considering the control procedure which allows the fixed point to change, it is clear that reducing the strength of the couplings proves beneficial for the system stabilization, in agreement with <cit.>. The same observation holds for the other classes of interactions (data not shown). As a final check, we show in Figure <ref> the new fixed points as obtained with the stabilized matrix A', compared to the trivial fixed point obtained when the matrix of couplings is switched off, i.e. assuming that each species is bound to evolve on a isolated niche.  The method here devised yields a genuinely complex stable equilibrium, which appears to be shaped by the couplings established among interacting units, despite the global tendency to reduce their associated strengths, as revealed in Figure <ref>.   
	
	
	
	
	
	
	
	
	
	
	

§ CONCLUSION

	
	
	In ecological systems, species interact in pairs and according to well-defined modalities, as e.g. predator-prey, mutualistic or competitive. 
	The dynamics of an ecosystem can be ideally modeled by resorting to a collection of ordinary differential equations for the evolution of the associated mean-field concentrations. This amounts to operating under the deterministic approximation, by deliberately neglecting the role played by finite size fluctuations and disregarding spatial variability. Each
	governing differential equation combines two distinct contributions: (i) local terms, also called reactions, forged to mimic the self-dynamics of the examined species and (ii) non local terms, which stem from the couplings among species. The matrix of interactions is a table where all existing pairwise interactions are stored and categorized. Given a specific ecosystem, the reaction parameters of the model and the associated network of inter-species couplings, it is interesting to speculate on the conditions that allow for the existence of a stable equilibrium. Stability relates to resilience, the ability of the system to withstand changes that would alter its dynamical equilibrium. Working in this context, we here contributed with a novel approach, mathematically grounded on first principles,  to help identifying the topological features that should be possessed by a generic ecological network so as to ensure stability,  hence resilience. To this end we consider a system made up of N interacting species and assume the reaction dynamics to be of the logistic type. We then generate a network of interactions, encoded in a weighted adjacency matrix, which is prone to instability. By rewiring the assigned links, and their associated weights, we can drive the system stable, via two alternative strategies which preserve, or not, the initial fixed point. The method consists of a generalized version of the topological control scheme designed in <cit.>, where diffusive (hence, linear) couplings were instead assumed to drive the exchanges between species. The technique here developed implements minimal modifications to the spectrum of the Jacobian matrix responsible for the stability of the underlying equilibrium and trace these changes back to species-species interactions. The analysis carried out within this operating framework, suggests that  predator-prey interactions exert a stabilizing effects, in qualitative agreement with the conclusion reached in <cit.>. Furthermore,  it has been found that a preponderance of weak interactions is beneficial to stability <cit.>.
	
	
	
	
	
	

§ COUPLINGS IN THE ORIGINAL SYSTEM

	
	
	
We shall here briefly discuss the algorithm that we have employed to generate the initial matrix of interactions A. For each pair i,j (with i and j running from 1 to N) we draw a random number p from a uniform distribution in the interval [0,1]. If p>p_0 a link that goes from j to i is established. Here 0<p_0<1 is a free parameter.  
The strength of the link, namely the element A_ij of matrix A,  is assigned as follows: we extract a random number from a uniform distribution defined in the compact interval [-c_0, 1-c_0], with c_0>0, and multiply the selected number by a scalar amplitude factor a. It is immediate to prove, that increasing c_0 makes the system progressively more unstable. This follows a direct application of the Gershgorin theorem, mentioned in the main body of the paper. The relative abundance of the pairs, typified as in the main text, can be on average computed and shown to be related to the choices of p_0 and c_0. The percentage of null entries (no links) of A will be 1-p_0, the percentage of negative entries p_0c_0, while the percentage of positive entries is p_0(1-c_0). Building on this observation, we obtain the following estimate for f_(·,·), the frequencies of occurrence of the different classes:	
	
    f_(0,0)   =    (1-p_0)^2
    
    	f_(0,-)   =    2p_0(1-p_0)c_0
    
    	f_(0,+)   =    2p_0(1-p_0)(1-c_0)
    
    	f_(+,-)   =    2p_0^2c_0(1-c_0)
    
    	f_(+,+)   =    p_0^2(1-c_0)^2
    
    	f_(-,-)   =    p_0^2c_0^2.

Each quantity is to be multiplied for a factor (N^2-N)/2 to obtain the abundances of elements, which are numerically tested in Figure  <ref>. The red crosses depicted in the Figures enclosed in the main body of the paper refer to equations (<ref>).


	
	

	
	
	

§ EXTENDING THE ANALYSIS TO ACCOUNT FOR A GENERALIZED NONLINEAR REACTION TERM.

	
	The control scheme developed in this paper can also be applied to a more general reaction model. Assume that the logistic dynamics is replaced by a generic nonlinear function f(x_i,r_i):
	
	
    ẋ_̇i̇=f(x_i,r_i) + x_i∑_j≠ iA_ijx_j

	where r_i identifies an arbitrary set of constant parameters.

	The associated fixed point is obtained by solving: 
	
    f(x_i,r_i) + x_i∑_j≠ iA_ijx_j=0.

	By performing the change of variables y_i≡ x_i/x_i^* and B_ij=A_ijx_j^* we obtain:
	
    ẏ_̇i̇=1/x_i^*f(y_ix_i^*,r_i) + y_i∑_j≠ iB_ijy_j

	
	In the new variables, the fixed point equation takes the form:
	
    1/x_i^*f(x_i^*,r_i) + ∑_j≠ iB_ij=0
 
	where use has been made of the condition y_i^*=1.

	The stability analysis requires introducing a modest perturbation v_i around the fixed point y_i^*=1, which amounts to writing y_i=1+v_i. The inhomogeneous perturbation v_i will evolve according to: 
	
    v̇_i = v_i[1/x_i^*∂ f/∂ y_i(y_ix_i^*,r_i)|_y_i=1 + ∑_j≠ iB_ij] + ∑_j≠ iB_ijv_j≡∑_j C_ijv_j.

	where the matrix C is defined as
	
    C_ij≡{[                           B_ij                      if i≠ j; ∑_k≠ iB_ik + 1/x_i^*∂ f/∂ y_i(y_ix_i^*,r_i)|_y_i=1  if i=j,;                                                             ].

	and the stability of the fixed point is ultimately controlled by the sign of the real part of the eigenvalues of C.
	In analogy with the procedure described in the main text, we then define D as
	
    D_ij≡{[           B_ij       if i≠ j; 1/x_i^*f(x_i^*,r_i)  if i=j,;                              ].

	which is a zero-row-sum matrix, because of the fixed point condition, ∑_jD_ij=∑_j≠ i B_ij+1/x_i^*f(x_i^*,r_i)=0.

	We can then modify the matrix D, which hence transforms into D', so as to enforce the desired stability, following the recipe outlined in the main body of the paper. The implemented changes can be interpreted as follows:
	
	
  * for the off-diagonal entries (i≠ j) we impose D'_ij = B'_ij = A'_ijx_j^*' which enables to calculate the elements of the controlled interaction matrix. As usual, we can decide to freeze the fixed point to its original value or modify it consistently, while assuming constant the parameters of the model.
		
  
  * for the diagonal entries (i=j) we impose D'_ii=1/x_i^*'f(x_i^*',r_i') where r_i' is the new vector of parameters of the stabilized system.  This latter can be readily obtained by inverting the above equation (f need therefore to  be invertible, with respect to r). Notice that, in general, only a subset of the elements of r_i, need to be adjusted. As remarked above, it is alternatively possible to leave the parameters r_i  unchanged, and modify x^*
 
	

	
	


§ THE CONTROLLED MATRIX D' IS REAL AND ZERO-ROW-SUM AS D IS

	
	Recall the definition of D':
	
    D'=D+Φ( δΛ)Φ^-1≡D+δ D

	where Φ is the matrix whose columns are the eigenvectors (ϕ^(1),...,ϕ^(N)) of D. Observe that since D is real and zero-row-sum, it is sufficient for our purposes to prove that δ D exhibits the same properties.

	
	The elements of D' are real.  Consider the generic entries of δ D:
	
    (δ D)_il=∑_jΦ_ijδΛ^(j)(Φ^-1)_jl.

	Isolate the real and imaginary parts of every element in Eq. (<ref>). It is immediate to see that, being δΛ^(j) real ∀ j, the imaginary part of (δ D)_il reads:
	
    (δ D_Im)_il = ∑_j δΛ_Re^(j)[(Φ_Re)_ij(Φ^-1_Im)_jl+(Φ_Im)_ij(Φ^-1_Re)_jl].

	
	To match condition (δ D_Im)_il=0, the term on the right hand side of Eq. (<ref>) should be zero. To prove this, let us first recall that the eigenvalues of a real matrix are either real or complex and come in conjugate pairs. Consider a complex eigenvalue Λ^(α), by definition  Dϕ^(α)= Λ^(α)ϕ^(α). Taking the complex conjugate yields 
	D (ϕ^(α))^*=(Λ^(α))^* (ϕ^(α))^* where (·)^* stands for the complex conjugate and where we have used the condition  D=  D^*. Let  the indexes α and β to be defined by the  relation (Λ^(α))^*=Λ^(β). We can immediately conclude that ϕ^(β)=(ϕ^(α))^*. Hence
	
	
    [  (Φ_Re)_iα=(Φ_Re)_iβ; (Φ_Im)_iα=-(Φ_Im)_iβ ]

	for each node index i and every pair (α,β) of complex conjugate eigenvalues. Consider now the equation
	
    (ϕ^-1)^(α)D=Λ^(α)(ϕ^-1)^(α)

	with (ϕ^-1)^(α) left eigenvector of D, corresponding to the α-th row of Φ^-1. Proceeding in analogy with the above, one eventually gets:
	
    [   (Φ^-1_Re)_α l=(Φ^-1_Re)_β l; (Φ^-1_Im)_α l=-(Φ^-1_Im)_β l. ]

	Return now to Eq. (<ref>). By performing the summation on j=α and j=β, using Eq. (<ref>) and Eq. (<ref>) and the fact that the corrections δΛ^(α) and δΛ^(β) are complex conjugated like the original eigenvalues  Λ^(α) and Λ^(β), we finally conclude that the terms of the sums in Eq. (<ref>) cancel in pairs relative to complex conjugate eigenvalues.

	The remaining terms in the summation correspond to real eigenvalues. Without loss of  generality, we can always consider the relative left and right eigenvectors to be real, up to a constant scaling factor. In detail, if Λ^(γ)∈ℝ, (Φ_Im)_iγ and (Φ^-1_Im)_γ l can always be set equal to zero ∀ i,l. This implies that all the remaining terms in equation (<ref>)  disappear, one by one. Then, summing up, (δ D_Im)_il=0 ∀ i,l.

	
	
	D' is a zero-row-sum matrix.
	 Matrix D zero-row-sum. Hence, the vector 1 (with all entries equal to one) is the right eigenvector of D corresponding to Λ^(1)=0:
	
    ∑_jD_ij=0      D1=0.

	Recall that the proposed approach, implies a modification of the eigenvalues of matrix D, while keeping the eigenvectors unchanged. As a consequence,  vector 1 is still solution of the eigenvalue problem. Moreover, the zero eigenvalue is not responsible for the instability and it is therefore preserved upon application of the control scheme. Hence:
	
    D'1=0    ∑_j(D')_ij=0

	which proves the claim.
	
	
	
		

§ ON THE CONDITIONS OF CONTROLLABILITY.

		
		
As explained in the main text, the proposed control method is based on shifting the eigenvalues of a zero-row-sum matrix in the complex plane by applying to their values a real and negative correction (which is identically equal to zero, for the subset of eigenvalues which should be preserved). The goal of this Appendix is to single out the conditions which allow for the control procedure to be effectively implemented.

		
Denote by ℳ the set of indices corresponding to the eigenvalues which are responsible for the instability and which should be modified by the controller.  Denote by  δΛ_i,  i∈ℳ, the (real) shift imposed to the selected ensemble of eigenvalues for stabilization. In the sequel, we will assume that the translation 
takes the interested eigenvalues to a constant  value, R, smaller than r_min (other strategies can clearly be adopted,  consequently altering the analysis reported below):

		
    (Λ_Re)_i+δΛ_i=R<r_min

		In choosing the constant R we must ensure that the applicability constraint (<ref>) is satisfied:
		
    s̃_i-∑_jΦ_ijδΛ_jΦ_ji^-1≥0  ∀ i

		Imposing condition (<ref>),  one gets:
		
    s̃_i - R∑_j∈ℳΦ_ijΦ_ji^-1 + ∑_j∈ℳΦ_ij(Λ_Re)_jΦ_ji^-1≡ k_i - R∑_j∈ℳΦ_ijΦ_ji^-1≥0 ∀ i.
 
		where we defined k_i=s̃_i+ ∑_j∈ℳΦ_ij(Λ_Re)_jΦ_ji^-1.

		
		
		
		If  k_i - r_min∑_j∈ℳΦ_ijΦ_ji^-1 is positive, for each index i ∈ [1,N], then it is sufficient to set R=r_min for achieving stabilization. The complementary situation, in which the above quantity turns out to be negative for some i, is more intricate, as we shall clarify hereafter. 
				
		First, it is convenient to sort the node indices in such a way that for indices i ∈ [1,n], with n<N, the quantity ∑_j∈ℳΦ_ijΦ_ji^-1 is positive, while it takes negative values for the remaining indices of the collection, i ∈ [n+1,N]. 
		Recalling that R is bound to be smaller than r_min, the following inequalities hold:
		
    k_i - R∑_j∈ℳΦ_ijΦ_ji^-1≤ k_i - r_min∑_j∈ℳΦ_ijΦ_ji^-1∀ i ∈ [n+1,N]
 
		
    k_i - R∑_j∈ℳΦ_ijΦ_ji^-1≥ k_i - r_min∑_j∈ℳΦ_ijΦ_ji^-1∀ i ∈ [1,n]
 
		where R is chosen so as to make the term on the left hand side positive. The inequality  (<ref>) is obviously violated if the expression on the right hand side is negative for at least one i in the interval [n+1,N]. In this case, the system cannot be controlled, using the recipe here discussed (which amounts, among the other specificities, to select a constant R).  For what concerns the other inequality (<ref>), the righthand term is instead allowed to be negative. Suppose that this happens for indices i ∈ [1,ñ] with ñ≤ n (the indices are imagined to be properly sorted). Then, the constant R must be smaller than r_min, let us say R = r_min - ϵ, with ϵ>0. Substituting it into (<ref>) we obtain a lower bound for ϵ:
		
    ϵ≥max_i∈[1,ñ](r_min - k_i/∑_j∈ℳΦ_ijΦ_ji^-1).

		
		Recall again that  R is (our choice) constant.  Assuming R = r_min - ϵ, one could eventually loose the controllability condition for indices i ∈ [n+1,N]. The following additional condition needs therefore to be considered:
		
    k_i - r_min∑_j∈ℳΦ_ijΦ_ji^-1 + ϵ_min∑_j∈ℳΦ_ijΦ_ji^-1≥ 0

		thus resulting in  an upper bound for ϵ: 
		
    ϵ≤min_i∈[n+1,N](r_min - k_i/∑_j∈ℳΦ_ijΦ_ji^-1).

		Another necessary condition for controllability is then found by imposing that the upper bound for ϵ is larger than the lower bound (both positive):
		
    0<max_i∈[1,ñ](r_min - k_i/∑_j∈ℳΦ_ijΦ_ji^-1) 
    		≤min_i∈[n+1,N](r_min - k_i/∑_j∈ℳΦ_ijΦ_ji^-1)

		
		which coincides with the equations reported in the main body of the paper.
		
	

	
	
	
	
	

§ CONTROLLING WITHOUT MODIFYING THE FIXED POINT

	
	 The Figures appearing in this Appendix illustrate the results obtained for the control strategy where the fixed point is kept unchanged while the inverse carrying capacity vector s is modified together with the weights of the interaction matrix. Figures <ref>,  <ref> and <ref> represent the analogue of Figures, respectively  <ref>, <ref> and <ref>, already discussed in the main body of the paper, which are obtained with the alternative control strategy.
		
	
		
	
	
	
	
	
	
	
	
	

	
	

tocchapterBibliography

ieeetr 




