A geometrical approach to control and controllability of nonlinear dynamical networks

In spite of the recent interest and advances in linear controllability of complex networks, controlling nonlinear network dynamics remains an outstanding problem. Here we develop an experimentally feasible control framework for nonlinear dynamical networks that exhibit multistability. The control objective is to apply parameter perturbation to drive the system from one attractor to another, assuming that the former is undesired and the latter is desired. To make our framework practically meaningful, we consider restricted parameter perturbation by imposing two constraints: it must be experimentally realizable and applied only temporarily. We introduce the concept of attractor network, which allows us to formulate a quantifiable controllability framework for nonlinear dynamical networks: a network is more controllable if the attractor network is more strongly connected. We test our control framework using examples from various models of experimental gene regulatory networks and demonstrate the beneficial role of noise in facilitating control.


Supplementary Tables
Supplementary Table 1: Edge control in the T-cell system. The first and the second columns, respectively, give the executor and the receiver of the coupling edges in the T-cell system. For control associated with each edge, the minimum control time τ m exhibits a power-law scaling behavior: τ m = α · |µ c − µ n | β , where µ c is the critical coupling strength (the third column). The fourth and the fifth columns list the scaling parameters α and β obtained by least squares fitting, with mean squared error (MSE) given in the sixth column. Note that Apoptosis represents the outcome of cellular signaling and it has inhibitory regulations to all the other nodes in the network. If Apoptosis is activated, the system will reach the desired normal state. Because of this the node Apoptosis appears in every row of the first column.    We also test the case of driving the system from a normal state to a cancerous state. In particular, each of the two experimentally adjustable parameters, the edge from node "Caspase" to "Apoptosis" and the self edge of "Apoptosis," can be perturbed for the control. We find that, for each control parameter, the relationship between its strength and the minimal control time also follows an algebraic scaling law with the scaling parameters: µ c ≈ 0.8778, α ≈ 1.1626, and β ≈ −0.6268 [equation (2) in the main text]. The control result is in accordance with the clinical studies revealing that the T-LGL leukemia disease results from dysregulation of apoptosis [1,2].
Supplementary Note 2: Algorithms to find attractors. Given a nonlinear dynamical network, the following procedure can be used to locate all the attractors.
1. For a given parameter set, we define the search space according to the maximum and minimum possible values of each state variable. For a GRN, the maximal value of each steady state is the activation rate divided by the degradation rate whereas its minimum value is zero. For example, for our two-node GRN model, we obtain that the maximal initial value for x 1 is (a 1 + b 1 ) · k −1 under the assumption that the leakage term is negligible.
2. Divide the phase space into a grid to generate a large number of initial conditions (grid points). Evolve each initial condition under the system dynamics to determine the final attractors of the system. Increase the grid resolution until no new attractors appear. For example, for our two-node GRN, the initial conditions are chosen from a 11 × 11 grid in the two-dimensional phase space region determined by the respective ranges of the dynamical variables. There are then 121 different initial conditions, which lead to four distinct attractors. Doubling the grid resolution results in no new attractors, enabling us to conclude that there are four distinct attractors.
3. For each attractor, calculate the eigenvalues of Jacobian matrix to determine its relative stability.
Supplementary Note 3: Parameter control method in stochastic model. Our nonlinear control framework is also applicable to stochastic systems. To demonstrate this, we convert the ODE model of two-node system [Eq.
(3) in the main text] to a set of Langevin equations for biochemical reactions and use the Gillespie algorithm to approximate the stochasticity [1]. Supplementary Fig. 1 shows a particular control process from attractor A to B. As discussed in the main text, A has low abundance in both x 1 and x 2 , and B has high abundance in x 1 and low abundance in x 2 . From the attractor network in Fig. 5(d) of the main text, we see that by increasing a 1 , we can drive the system from A to B. Thus, we first set the system in A with a 1 = 1.0. At t = 60, we increase a 1 to 1.5, and at t = 120, we change a 1 back to 1.0. We see that, for t > 120 when the perturbation on a 1 has been withdrawn, the system spontaneously evolves into the attractor B.
Supplementary Note 4: Control based on multiple parameters. To figure out the optimal parameter combinations for controlling nonlinear networks of even moderate size is computationally prohibitive at the present. However, for small networks, this can be done. For example, for the three-node GRN system studied in this paper, multiple parameter control is needed to induce certain state transition, e.g., from attractor H to attractor B in Fig. 7 in the main text. However, for the T-cell network, there are several dozens tunable edges. A surprising finding is that, because of the simplicity of the attractor network, it is only necessary to apply perturbation to any one of the 48 tunable edges to realize control. It is interesting to note that, the sequential combination of control dosage guided by the attractor network may be regarded as a kind of control based on parameter combination. For example, it has been known that the human embryonic stem cells (hESCs) can be differentiated to a pancreatic fate under stepwise exposures to different signaling factors [2,3]. However, the functional pancreatic fate would not occur if the differentiation steps are permuted. The non-interchangeable sequence can be understood using the concept of the attractor network, where one parameter modification represents one stage protocol. For example, in Fig. 5(d) in the main text, the hESCs correspond to attractor A. By increasing the value of a 1 (stage 1) and then that of b 2 (stage 2), we can drive the system to attractor C. However, if we increase b 2 first (stage 2 first), we will not be able to drive the system to attractor C by increasing the value a 1 (stage 1) [see also the top branch of attractor D in Fig. 5(a) in the main text]. We emphasize that these results do not imply that the sequential combination method is restricted to perturbation to one parameter at a time. Indeed, recent works suggested that simultaneous therapy with two drugs can be much more effective than sequential therapy [4].
Supplementary Note 5: Algorithmic complexity and computational cost. In the attractor network, a link is determined when the undesired attractor vanishes upon application of control parameter perturbation of certain magnitude (within some predefined range). For each attractor, we perform a bisecting search for all the tunable parameters to establish the possible links in the attractor network. In each search, the number of checks to see if the attractor has disappeared is log 2 (1/∆), where ∆ is the accuracy in the estimate of the critical perturbation amplitude. The total number of bisecting searches is the number of attractors, N S , multiplying the number of tunable parameters, N L . In a nonlinear system, the number of attractors depends on the system size. For example, for a boolean network, there are 2 N possible states. While the actual number of attractors can be much less than 2 N , it depends on the network size N . The number of control parameters is the number of controllable links. Assuming that the network is sparse (as in many biological networks), the total number of checks is ≈ log 2 (1/∆) × N S × ρN 2 , where ρ is the connection probability. The bisecting searching needs to be performed once for a sufficiently large time duration τ of the control. The reason is that, if the control is achievable for a link for a longer duration, the same control can be realized for a shorter duration but with a larger control amplitude.
For nonlinear and complex dynamical systems, the relationship between control perturbation and escaping from an attractor can be discontinuous and/or non-monotonic. In such a case, a "blind" application of the bisecting search to build up the attractor network and to estimate the computational cost may not be effective or even fail. This difficulty, however, can be overcome by using the method of parameter continuation with multiple initial conditions, which is standard in generating and analyzing the bifurcation diagrams of nonlinear dynamical systems. Especially, one can choose a small number of random initial conditions and a small set of parameters in a physically/biologically meaningful range, and determine the distinct attractors that the system possesses. The result will be a global picture of the possible attractors of the system in a small number of parameter intervals, which will facilitate greatly the computational task.
Supplementary Note 6: Phase diagram. A phase diagram illustrates how different choices of the parameters affect the system's stable asymptotic states (attractors). Supplementary Fig. 2 shows the phase diagram of the two-node GRN system with respect to variations in a 1 and b 2 . Each point in the diagram represents, for the particular combination of values of a 1 and b 2 , the whole set of possible attractors of the system. There are seven possible combinations of attractors in the diagram. As shown in Fig. 5(d) in the main text, if we set out to control the system from attractor A to attractor C, we can first steer the system from A to B through perturbation on a 1 and then drive the system from B to C by tuning b 2 . Note that the sequence of a 1 and b 2 is important here: if we first perturb b 2 , the system in attractor A will be driven to attractor D. When this happens, the system will remain in D, regardless of any additional parameter adjustments, rendering unrealizable the control goal to drive the system into attractor C. This scenario can also be seen from the phase diagram in Supplementary Fig. 2. In particular, starting from the dark red region, with the particular combination of parameters, there are 4 attractors: A, B, C and D. As we increase a 1 and set the parameters to the red area, attractor A disappears and its basin is merged into that of one of the remaining three attractors (B, C, and D). From Fig. 4(e) in the main text, we see that A migrates into B. Thus, when we withdraw the perturbation on a 1 and turn on the perturbation on b 2 , the possible attractors of system are C and D (the light green area). Figure 5(d) in the main text also indicates that B migrates to C. We note that, in the phase diagram, it is difficult to distinguish into which state A merges, but this can be readily accomplished through the attractor network in Fig. 5(d) in the main text, which indicates unequivocally that the basin of A is absorbed into that of D.
Supplementary Note 7: Termination criteria and calculation of control strength and minimum control time. For all systems studied in this paper, of particular importance to control is the relationship between control strength and minimum control time. The attractor network needs to be constructed first based on the procedures described in Methods in the main text. From the attractor network, we can obtain the control parameter for the transition between any two states. For each control parameter, we choose several different values which can realize the control and use a bisecting search to find the minimum control time. We then apply linear regression to the data on a double logarithmic scale to find the value of µ c . In the bisecting search process, the time interval is chosen to be t s = 10 −2 . For the T-cell and two-node GRN systems detailed in the main text, we choose T f to be 100. For the three-node GRN system, we choose T f to be 1000. The results for each system are illustrated in Supplementary Table 1 Table 3 illustrates the critical coupling strength and the minimal control time to realize control.