Co-evolutionary control of a class of coupled mixed-feedback systems Open Access

This work tackles the challenge of stabilizing a desired oscillatory pattern in a network with a fixed structure. We propose two distributed controllers to achieve this. These controllers are fundamentally different: the first is highly robust and eliminates nonlinearities but requires full knowledge of the reference system. The second, inspired by the neuromodulation of biological systems, only requires local error information. Moreover, the controllers are implemented as co-evolutionary rules for edges connecting an oscillatory node with the controlled system. This effectively results in an adaptive network where the adaptation rules shape the dynamics of the controlled nodes. We showcase the co-evolutionary controllers by exploring the synchronization of the controlled system, showing that our approach is independent of the network structure, and testing the performance of the controllers in time-varying networks. Our theoretical and numerical results show two distributed strategies to induce a desired oscillatory pattern by adaptively modulating node> dynamics.

Inducing and regulating a desired oscillatory behavior in networked systems is a fundamental problem in control theory, with ample relevance for both natural and engineered systems. Oscillations are ubiquitous in many real-world phenomena including industrial applications,^1–3 chemistry,^4,5 neural networks,^6,7 psychology,^8,9 and biology.^10,11 More specifically, in biological systems, rhythmic signals govern important physiological processes, including the heartbeat’s regulation, circadian rhythms, central pattern generators, and neuronal signaling.^12–15 From the engineering side, several technologies also need such periodic signals, for example, to achieve the coordination of automation processes, including robotic locomotion,^16,17 or in the developments of synthetic biological systems.^18–20 

Several naturally occurring oscillatory patterns are self-regulated and emerge from the intrinsic dynamics of the system. In contrast, there are many scenarios, both natural and engineered, where external control is required to induce or modify oscillations in a desired way. This is, for instance, relevant in neuromorphic engineering,²¹ where particular oscillatory patterns of biologically inspired neuronal networks are critical, for example, to mimic cognitive functions, coordinate sensorimotor functions, synchronize spiking neural networks, and enable efficient signal encoding, to name a few.^22–24 While our interests in this paper are mainly theoretical, designing controllers that can reliably induce desired oscillations could lead to advances in neuromorphic engineering and design, brain-inspired control systems, and adaptive network dynamics, among others.

In this paper, we design two controllers that impose a desired oscillatory pattern on a networked system without modifying its internal dynamics or interconnections. Our focus is on a class of coupled mixed-feedback systems where the interplay of (fast) positive and (slow) negative feedback plays a key role in the generation of sustained oscillations [see (1) and further details below].

The manuscript is organized as follows: in Sec. II, we motivate and describe the problem at hand and propose two kinds of controllers. The first controller, while robust, needs complete information from the reference. In contrast, the second one only requires local error information. Next, in Sec. III, we embed the designed controllers into a networked framework, implementing them as a co-evolutionary, or adaptive, rule of the edges connecting a node (the controller’s node) to the nodes of the controlled system. We compare and evaluate the controllers’ performance in three key scenarios: synchronization, arbitrary network structures, and time-varying network structures. We conclude in Sec. IV, where we also discuss potential extensions of our work.

Notation: $R$ and $C$ denote the fields of real and complex numbers, respectively. We use $sgn$ for the sign function and $O$ to denote the big-Oh order symbol. An adjacency matrix is denoted by $A$⁠, and a superscript $r$ (for “reference”) or $p$ (for “plant”), as in $A r$ and $A p$⁠, is used whenever it is necessary to distinguish the adjacency matrices of the reference and of the plant. Further notation is clarified when necessary.

As in (2), a key feature of mixed-feedback systems is the interplay between positive and negative feedback loops across different timescales, which is commonly observed in many biological models, particularly in neuronal dynamics. For example, it is argued²⁶ that mixed feedback is the fundamental mechanism for excitability and explains how biological systems are able to exhibit sustained and robust oscillations.^27,28 Moreover, mixed-feedback architectures have been recently exploited to develop control methods for neuromorphic systems.²⁵ 

Our main motivation to consider systems given by (1) is the theory developed in Ref. 29. Relevant to us is that they propose a procedure that allows one, for example, to design an adjacency matrix such that the output $x=( x 1,…, x n)∈ R n$ corresponds to a rhythmic profile of the network (defined as the inverse problem in Ref. 29). While brief, let us be more precise: a rhythmic profile is an $n$-tuple $( σ 1 e ı ϕ 1,…, σ 1 e ı ϕ 1)∈ C n$⁠, where $σ i∈ R$ represents the amplitudes and $ϕ i∈[0,2π)$ phases. If for all $i=1,…,n$⁠, the solutions $x i(t)$ of (1) converge to some periodic function with amplitude $σ i$ and phase $ϕ i$ defined as oscillating function,²⁹ then one says that the network is rhythmic. Assume further that $σ 1>0$ and that $σ 1≥ σ i$⁠, for all $i=2,…,n$⁠. The $n$-tuple $(1, ρ 2 e ı θ 2,…, ρ n e ı θ n)$⁠, where $ρ i= σ i σ 1$ and $θ i= ϕ i− ϕ 1$ for all $i=2,…,n$⁠, is called a relative rhythmic profile. By following the procedure summarized below, one can construct an adjacency matrix such that the solutions of (1) (locally) converge to a particular rhythmic profile:

1. Let $ω x=(1, ρ 2 e ı θ 2,…, ρ n e ı θ n)$ with $θ i≠ { 0 , π }mod2π$⁠.

2. Pick $μ 1=a+ıb$ with $a,b>0$⁠. Further choose $μ 3,…, μ n∈ R$ with $μ i<a$ for all $i=3,…,n$⁠. Define $D=diag⁡( μ 1, μ 1 ¯, μ 3,…, μ n)$⁠.

3. Find an invertible matrix $Q= [ w x w x ¯ B ]∈ C n × n$⁠, where $B∈ R n × ( n − 2 )$⁠.

4. Define $A=QD Q − 1$⁠.

We now consider the following problem: suppose we are given a reference rhythmic network and a plant network, both of type (1). The adjacency matrix of the plant cannot be adjusted, but we want to render the plant rhythmic with respect to the reference’s profile. Since the adjacency matrix of the plant is fixed, we shall design a controller that solves our problem.

In the following, we propose two different kinds of controllers that render the origin of (4) locally stable. The first one eliminates the sigmoidal nonlinearities, which results in a robust controller but requires complete knowledge of the reference. For comparison purposes, we propose an additional controller inspired by the synaptic modulation of neuronal systems. Such a controller does not require any knowledge of the reference. For clarity of our exposition, these two types of controllers are first provided without considering their implementation, see Propositions 1 and 2. After that, we propose a particular way of implementing them in a networked system where all nodes are of type (1), including the controller. This will mean that the controller action is implemented via the adaptation of the weights connecting the controller node and the plant’s nodes; see more details in Sec. III.

We describe some generalities about the simulations that we present below. When required, more details are provided, and all codes are available in Ref. 30. In our simulations, the function $S$ is chosen as the odd sigmoid $S(⋅)=tanh⁡(⋅)$⁠, the small parameter $ε$ is fixed at $ε= 1 100$⁠, while many of the other parameters are chosen at random. When we say that a parameter is chosen “at random,” we always mean it under a uniform distribution and specify the interval of possible values. The following steps align with the algorithm described inSec. I:

1. The relative amplitudes $ρ i$⁠, $i=2,…,n$⁠, are chosen at random within the interval $( 1 3,1)$⁠. The lower value $1 3$ is simply chosen so that the $i$th relative amplitude is not too small in the simulations. The phases $θ i$⁠, $i=2,…,n$⁠, are chosen at random within the interval $(0,2π)$⁠. The probability of $θ i=π$ is zero.

2. For the leading eigenvalue $μ 1=a+ıb$⁠, we let $a=1$ and $b$ is randomly chosen within the interval $( 1 100, 1 10)$⁠. The choice of $b$ is so that the period of the rhythmic profile, which is²⁹  $T≈2π ( β ℑ ( μ 1 ) ) − 1$⁠, is not too small. The rest of the eigenvalues $μ i∈ R$⁠, $i=3,…,n$⁠, are chosen at random within the interval $(0, 9 10a)$⁠. This provides $D=diag⁡( μ 1, μ 1 ¯, μ 3,…, μ n)$⁠.

3. The complement matrix $B$ is a randomly generated sparse matrix with nonzero coefficients within the interval $(0,1)$⁠. This matrix is generated until $Q= [ w x w x ¯ B ]$ is invertible.

4. The adjacency matrices $A r$ and $A p$ are, thus, given via the previous algorithm as $A ∙=QD Q − 1$⁠. Since the algorithm to obtain the adjacency matrices is mostly random, both matrices are, with probability $1$⁠, different. Likewise, since $B$ is sparse and randomly generated, the topologies of the adjacency matrices are, with probability $1$⁠, different. We emphasize that, from the results we present below, the plant’s adjacency matrix $A p$ does not have to be rhythmic. To keep all systems (reference, plant, and controller) within the same context, we have opted to keep it rhythmic for the simulations.

In addition, and due to the choice of $μ 1$⁠, $β ∙$ is randomly chosen within the interval $(0,1)$⁠, and the corresponding $α ∙$ is set as $α ∙=1+ε− β ∙+ 1 100$ (see Ref. 29, Theorem 1).

To better see the effect of the controllers, we present our simulations in the following way: for the first $300$ time units, the controller is off, and so we see the open-loop response. At $t=300$⁠, the controller is turned on, and so from thereon, we see the closed-loop response. In our algorithms, this is realized by multiplying the controller by $l(t− t c)$⁠, where $l(x)= 1 1 + e − x$ and $t c=300$ is the time at which the controller is turned on. Since we chose some parameters at random, the plots we show are representative of several simulation runs.

A simulation of (8) with $F i$ given by (7) is provided in Fig. 2; recall the general considerations described in Sec. II A. Figure 2 shows in the first two rows and Fig. 2(a) the time series for $x i(t)$⁠, $X i(t)$⁠, $u i(t)$ and the mean norm of the error $1 n | x ~ |$ in logarithmic scale for $n=4$ nodes. In this and all following simulations, we prefer to show $1 n | x ~ |$ because we keep all the controller constants the same, even for different numbers of nodes. Therefore, for comparison purposes, it is convenient to normalize the error. In Fig. 2(b), we show the mean error for a simulation similar to the previous one but for $n=100$⁠. Figure 2(c) also shows the mean error for a simulation with $n=100$ nodes but for time-varying adjacency matrices. More precisely, for this simulation, each entry $A i j ∙$ (for both the plant and the reference) of the adjacency matrices is of the form $A i j ∙= A i j ∙(t)= A ¯ i j ∙(1+ 1 5sin⁡( ω i jt))$ where $ω i j$ is some random frequency within the interval $(0,1)$ and $A ¯$ is a random rhythmic matrix. Since the controller fully uses these adjacency matrices, its performance is also good. Finally, on Fig. 2(d), we show another simulation for $n=100$ nodes but now with a mismatch between the parameters used by the controller and those of the reference. More precisely, for this simulation, we perturb the parameters $α r$⁠, $β r$⁠, and $A i j r$ used by the controller in (6), by some small random number within the interval $(− 1 20, 1 20)$⁠. Such small perturbation is different for each parameter and each entry of $A r$⁠. We notice that since the mismatch is relatively small, the performance of the controller is still reasonable, although the error is roughly one order of magnitude larger than the error without the mismatch [Fig. 2(b)]. For all these simulations, we have set $k i=k=100$⁠. It follows from Proposition 1 that the larger $k$⁠, the smaller the error.

As already mentioned, a drawback of the controller presented in Sec. II B is its dependence on the full knowledge of the reference’s parameters. To offer an alternative, we propose a second controller whose behavior is inspired by the neuromodulation of synaptic weights in neuroscience and neuromorphic systems.

Intuitively, the controller $u i$ adapts to the rescaled error $k λ x ~ i$ (quickly if $τ$ is large) and effectively provides a proportional feedback $u i≈ k λ x ~ i$⁠. If the error $x ~ i$ is large, then the controller compensates the second sigmoid function to balance them out. If the error is small, then both sigmoid functions have roughly the same value and cancel out, which leads to the exponentially stable linear error dynamics (5). Besides the well-known advantages of introducing the controller in integral form, the idea just described has its inspiration in some neurological mechanisms, where the dynamic equation for $u i$ is in “leaky form,” similar to the way biological neurons process information. Likewise, the introduced controller can be regarded as an adaptive synaptic weight that modulates the influence of the error. All these interpretations become more evident in Sec. III.

The previous idea is formalized in Proposition 2.

Figure 3 shows a simulation of (3) with the controller given by (10); recall the general considerations described in Sec. II A. The figure shows in the first two rows and Fig. 3(a) the time series for $x i(t)$⁠, $X i(t)$⁠, $u i(t)$ and the mean norm of the error $1 n | x ~ |$ in logarithmic scale for $n=4$ nodes. In Fig. 3(b), we show the mean error for a simulation similar to the previous one but for $n=100$⁠. Figure 3(c) also shows the mean error for a simulation with $n=100$ nodes but for time-varying adjacency matrices in the same way as those for Fig. 2. Since this controller does not use the parameters of the plant, there is no equivalent simulation to that of the Fig. 3(d) in Fig. 2. For all the simulations in Fig. 3, $τ=1$⁠, $λ=1$⁠, and $k=100$⁠. It follows from Proposition 2 that the larger the $k$⁠, the smaller the error.

Here, the biological inspiration is more evident: the adaptation rule (14) has two main components; one is a “pre-synaptic” activity-dependent term $x ~ i x c$⁠, which is regulated by the “post-synaptic” error and is reminiscent of a Hebbian-like adaptation (one could say that the learning rate is mediated by how far the plant’s node is from the reference’s), and an intrinsic decay $−λ a i$ ensuring that the values of $a i$ remain bounded.

In Subsecs. III A–III C, we present the simulations of (11), and we compare through different application scenarios, the adaptation rules given by (12) and (13) and (14). In addition to the general considerations for our numerical simulations described in Sec. II A, we add the following: in the time interval $t∈(500,550)$⁠, we add a constant random disturbance, taking values in the interval $(−1,1)$⁠, to each node, different for each node and independent of the node. We do this to numerically test how the controlled plant recovers from a perturbed state.

To avoid repetitions, we mention here once that in our numerical simulations of (11) with the adaptive implementation given by (12) and (13) (specifically, Figs. 4, 6, and 8), we have used $δ= 1 100$ and $τ c=100$⁠. Regarding the adaptive implementation given by (14) (specifically Figs. 5, 7, and 9), we have used $τ c=100$⁠, $k=100$ and $λ=1$⁠. We have kept these values the same across the different simulations for comparison purposes; of course, the performance of the adaptation rules can be improved if one tunes such parameters adequately.

Figure 4 shows a simulation of (11) with the adaptation rule given by (12) and (13). The observed spikes in the weight $a i$ are due to the discontinuous form of the adaptation and are modulated by $δ$ in (13). In contrast, Fig. 5 shows a simulation of (11) with the adaptation rule given by (14). For this case, we do not see anymore the aforementioned spiking behavior since $a i$ is smooth. However, we observe that the error tends to increase periodically, and this coincides with $x c$ crossing zero; as for $x c=0$⁠, the system is in open-loop. This issue is barely noticeable in the time series, though. We moreover notice that on both adaptation rules, the closed-loop systems recover well enough after the added disturbance in $t∈[500,550]$⁠.

To highlight the fact that our approach does not depend on the particular topologies of the reference’s and plant’s adjacency matrices, in this section, we present simulations with explicitly different adjacency matrices. Figure 6 shows a simulation of (12) and (13), where for $n=4$⁠, we have included the graphs of the reference and of the plant, which are different. As in the simulations of Sec. III A, the spikes in the time series of $a i$ are due to the discontinuous implementation (12) and are modulated by $δ$ in (12). For comparison, Fig. 7 shows a simulation of (11). Both simulations follow the same overall format we have used so far: for $t∈[0,300]$⁠, the controller is off, and hence, we see the open-loop response; at $t=300$⁠, the controller is turned on, and we see how the plant follows the reference, while for $t∈[500,550]$⁠, a constant random disturbance is added to each node. We notice from the time series that the performance of both controllers is reasonably good.

As a final application and comparison setup, we consider here that the adjacency matrices of both the reference and the plant are time-varying. We do this similar to what was presented in Sec. II. More precisely, for the simulations of this section, each entry $A i j ∙$ (for both the plant and the reference) of the adjacency matrices is of the form $A i j ∙= A i j ∙(t)= A ¯ i j ∙(1+ 1 5sin⁡( ω i jt))$⁠, where $ω i j$ is some random frequency within the interval $(0,1)$ and $A ¯ ∙$ is a random rhythmic matrix, meaning that the underlying topologies of each network are also different. While we do not perform formal analysis under this setup, both adaptation approaches seem suitable to handle this scenario. On the one hand, we have that adaptation rules (11) and (12) use the full knowledge of the adjacency matrices. On the other hand, the adaptive rule (11) accounts for the time-varying adjacency matrices from the fact that $a i$ is regulated by the error. Figure 8 corresponds to (10) with (11) and (12), while Fig. 9 uses (14). The description of the results is similar to Secs. III A and III B. In both cases, we can observe that, despite the adjacency matrices being time-varying, both adaptation rules perform well.

This paper has explored co-evolutionary control approaches to render a desired oscillatory pattern stable for a particular class of dynamic networks, keeping their internal properties fixed. We have developed two kinds of controllers: one that uses the full information from the reference and another that only requires local error information. From our analysis, both controllers seem to perform reasonably well. Their main difference is their potential implementation. While the first controller (see Proposition 1) performs observably well, it is hardly implementable not only because it requires full information about the plant but also because of its particular form [see (6)]. In contrast, the second controller (see Proposition 2) requires only local error information. We have also implemented such controllers as co-evolutionary rules in an extended dynamic network and tested their performance in three relevant scenarios. In all cases, again, the controllers perform reasonably well.

An immediate open problem to consider in the future is to extend our approach to more general coupled systems. As a concrete example, one may want to explore ideas similar to those developed here, but for coupled spiking neurons, whose models also have mixed-feedback loops. A challenge for these models is that the nonlinearity is not bounded, in contrast to (2). Another possibility is to consider that the reference and the plant are not of the same type. Hence, even if a system cannot inherently produce and sustain specific oscillations, a controller may impose them.

The work of L.G.V.P. was founded by the Data Science and Systems Complexity (DSSC) Centre of the University of Groningen, as part of a PhD scholarship. H.J.K. would like to acknowledge the financial support of the CogniGron research center and the Ubbo Emmius Funds (University of Groningen). The authors thank the anonymous reviewers for their thorough comments that helped to improve the manuscript.

The authors have no conflicts to disclose.

Luis Guillermo Venegas-Pineda: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Hildeberto Jardón-Kojakhmetov: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Writing – original draft (equal). Ming Cao: Conceptualization (equal); Formal analysis (equal); Project administration (equal); Supervision (equal).

The data that support the findings of this study are openly available in CoevolutionaryControl-2025, at https://github.com/hjardonk/CoevolutionaryControl-2025, Ref. 30.