Tractable nonlinear memory functions as a tool to capture and explain dynamical behaviors

Mathematical approaches from dynamical systems theory are used in a range of fields. This includes biology where they are used to describe processes such as protein-protein interaction and gene regulatory networks. As such networks increase in size and complexity, detailed dynamical models become cumbersome, making them difficult to explore and decipher. This necessitates the application of simplifying and coarse graining techniques to derive explanatory insight. Here we demonstrate that Zwanzig-Mori projection methods can be used to arbitrarily reduce the dimensionality of dynamical networks while retaining their dynamical properties. We show that a systematic expansion around the quasi-steady-state approximation allows an explicit solution for memory functions without prior knowledge of the dynamics. The approach not only preserves the same steady states but also replicates the transients of the original system. The method correctly predicts the dynamics of multistable systems as well as networks producing sustained and damped oscillations. Applying the approach to a gene regulatory network from the vertebrate neural tube, a well-characterized developmental transcriptional network, identifies features of the regulatory network responsible for its characteristic transient behavior. Taken together, our analysis shows that this method is broadly applicable to multistable dynamical systems and offers a powerful and efficient approach for understanding their behavior.

expansion Eq. (13) to find for the first term on the right hand side (RHS) of Eq. (5), where the two last terms arise by differentiating the product (x b − x * b ) f bs (x s , τ ) w.r.t. x s and we have not written the x s dependence of x b * for brevity. The second term on the RHS of Eq. (5) is the expectation of this, obtained by replacing x b by x b * , Putting the two together gives for the time evolution Eq. (5) of F s : For consistency with Eq. (13), we now linearize the square brackets again in x b − x b * . In the second term, we similarly replace R s (x s , x b ) by R s (x s , x b * ) as the remaining factor in this term is already linear. Comparing then with the time derivative of the original linearized formula Eq. (13) gives, after appropriate relabelling of indices, the equation for the evolution of f bs in time Eq. (15). From the above expansion, one sees that the matrix l bb in Eq. (15) takes the form The form Eq. (16) in the main text is obtained by using the identity The latter can be obtained by differentiating Eq. (9) with respect to x s . To obtain the actual memory function from Eq. (19) is now straightforward as we have already worked out the required expectation in Eq. (A2). The first term on the r.h.s. of Eq. (A2), which we had previously kept to make the ensuing linearization easier to see, actually vanishes because of Eq. (9), yielding Using again the identity Eq. (A5), we obtain our main result Eq. (20).

APPENDIX B: SOLUTION FOR F
In this Appendix, we find the solution f bs (x s , τ ) to the differential Eq. (15). We start by restating the latter as where we have used that the factor R s in Eq. (15) is just the effective drift v s defined in Eq. (11). As the equation is linear in f bs , its derivatives it can be solved using the method of characteristics (see, e.g., Ref. [41]). Calling the curve parameter for a characteristic u, the characteristic equations can be read off from Eq. (B1) as Setting an arbitrary integration constant to zero, the first of these gives τ = u. To solve Eq. (B3), we call φ v the flow generated by v(x s ), which is defined as the solution of the differential equations The solution of Eq. (B3) is then where x s 0 is the value at the beginning of the characteristic curve (u = 0); the minus sign in the second argument of φ v reflects the backward in time propagation in Eq. (B3). We note for later that, as a consequence of Eq. (B6), the solution values at u 1 and u 2 are related by Finally, the solution of Eq. (B4) is using the initial condition Eq. (14) at τ = u = 0. From Eq. (B4), we see that the matrix exponential appearing here must be time ordered, with earlier "times" u appearing to the right of later ones. It now remains to express f bs (u) in terms of x s (u) and τ (u) = u. We fix aû =τ and callx s = x s (û). Using Eq. (B7)  with u 2 = u and u 1 = u then shows that the x s -solution Eq. (B6) can be expressed in terms ofx s as and, in particular, Changing the integration variable to τ =τ − u and dropping the hats then gives the solution Eq. (18) announced in the main text. Note that as τ =τ − u , the time ordering of the matrix exponential, is such that the earlier τ are now on the left. The appropriate time-ordered matrix exponential is defined formally via its Taylor series, with 1 the identity matrix and the integration in the other terms running over the range 0 < τ 1 < . . . < τ n < τ.

APPENDIX C: MAPPING OF SELF-CONSISTENT MEMORY TO DIFFERENTIAL EQUATIONS
We show in this Appendix how to map the subnetwork equations with self-consistent memory, to a set of differential equations. The self-consistent memory termM s (t ) is given by Eq. (24), so can be written as It is then straightforward to check that where the second term arises from the t dependence of the matrix exponential. The m b (t ) can therefore be obtained numerically by integrating the differential Eq. (C5) together with the subnetwork equations with (self-consistent) memory: The appropriate initial conditions for the auxiliary variables follow from Eq. (C4) as m b (0) = 0.

APPENDIX D: CHANNEL DECOMPOSITION
We begin by writing the expression for the memory function explicitly, combining Eqs. (14), (18), (20), and (21): where the first three factors are evaluated at x s . We now swap index labels and group the sums into a more intuitive form: As discussed in the main text, the expression up to before the exponential represents a change in the deviation of the bulk species concentration x b from its QSS values over some small time interval, in response to changes in the subnetwork concentrations x s [see also Eq. (I4) below]. In the factor ∂R b /∂x s , only those bulk species b contribute whose time evolution depends explicitly on the subnetwork species s driving the bulk time evolution via R s . The b can then be interpreted as outgoing channels for the signal from s . After propagation in the bulk network, the signal returns via another bulk species. Here only bulk species b contribute that appear explicitly in the time evolution of subnetwork species s as indicated by the factors ∂R s /∂x b . The b can therefore be interpreted as incoming channels. Overall, we have memory effects from s onto s, via an outgoing channel (s to b ) and an incoming channel (b to s). Consistent with this interpretation, the outgoing channel "susceptibilities" ∂R b /∂x s are evaluated for the past, i.e., sending, state x s ≡ x s (t ) of the subnetwork. The incoming channel susceptibilities ∂R s /∂x b , on the other hand, are evaluated at the current time t as shown by the propagation via φ v across the time difference τ = t − t . Within the self-consistent approximation Eq. (C5), this  propagation corresponds directly to evaluation at the current state x s (t ).

APPENDIX E: SELF-CONSISTENT CHANNEL DECOMPOSITION
The channel decomposition of Sec. D can also be applied to the self-consistent memory approximation, as we now outline. Writing out the self-consistent memory term Eq. (C2) explicitly and reordering and relabelling terms as in Eq. (D4) gives where From this last representation, it follows that the m sbb s (t ) vanish at t = 0 and obey the differential equations: The channel-decomposed memory can therefore also be calculated from differential equations (see Appendix G). Of course, one only needs to find the m sbb s for combinations (sb) and (b s ), where the corresponding channel susceptibilities are nonzero.

APPENDIX F: EXACTNESS OF MEMORY
We show that the self-consistent memory m b (t ) is exact when both R s and R b contain at most linear terms in x b . In such a case, the full system can be written as x s ) and the QSS value x b * (x s ) is an arbitrary function of x s . We now want to show that thex b correspond exactly to the m b from the self-consistent ZMs method. To do this, we work out their evolution in time: By using that for an x b -linear system as assumed here, one has f 0 bs = ∂R s /∂x b , the above can be rewritten as Using then the definitions Eqs. (16) and (21), we obtain an expression equivalent to Eq. (C5), thus showing thatx b = m b when we start from the same initial conditionx b = 0, i.e., the bulk at QSS. We test the above exactness statement on two different examples that have a linear dependence on a particular species but nonlinear dependencies on other species: a minimal bistable system as described in Ref. [30], and the Brusselator, which is capable of achieving limit cycles [31]. As expected from the above derivation, the self-consistent memory captures the behavior of both systems exactly (Figs. 9 and 10). As further shown in Fig 10, the original nonlinear projection method ZMn is also accurate at capturing the dynamics though not necessarily exact. The corresponding memory function is shown in Fig.11. (We note that the memory functions of the Brusselator grow exponentially in a way that forces memory terms to cancel out to zero at the fixed point; this leads to numerical challenges that we do not pursue here.)  16) and (14) that l bb and f 0 bs are both constant, i.e., independent of x s . Accordingly [compare Eqs. (18),(20), and (24)]. the memory functions of the ZMn and ZMs projections also become identical, and the corresponding channel decompositions are also the same.
To illustrate the linearized dynamics approach, we perform a channel decomposition of the amplitude (value at τ = 0) of the linearized memory in the neural tube system as we did in Ref. [27], but now for the method derived in this study (Fig. 12). We find similar profiles to those found in Ref. [27]. The results highlight the relative weakness of the memory from Olig2 into Nkx2.2 via Pax6, supporting the conclusions of Ref. [27]. In addition, channel decomposition for the ZM trajectory shown in Fig. 6 provided memory terms affecting Nkx2.2 and Olig2 (Fig. 13). The method derived in this study is, however, significantly more powerful as it does not rely on an expansion near a steady state and gives access to the full memory and its channel decomposition as described in Appendixes D and E.

APPENDIX H: COMPARISON WITH ALTERNATIVE MEMORY FUNCTION APPROXIMATION
Gouasmi et al. [18] propose an approximation for the memory function for the case where the projection Eq. (10) is defined not by setting the bulk coordinates to their x sdependent QSS values, but simply to zero: FIG. 11. Memory functions for the system detailed in Ref.
[30] with the parameters chosen for Fig. 10. (a) Using the ZMn. (b) Using the method from Gouasmi et al. as extended to QSS projection in Appendix I. The x axis shows the concentration of the subnetwork species x 1 while the y axis indicates time difference τ . By construction, the two memory function approximations predict the same value (scale bar to the right) at τ = 0, as they only differ in how they propagate the memory over time. At τ > 0, the memory functions are relatively similar for x 1 ∈ [0, 6] but become progressively different as x 1 grows beyond this range; for x 1 10, the G-QSS method predicts a negative memory function for all τ that leads to its poor performance as observed in Fig. 10. The function F s (x s , x b , τ ) still evolves according to Eq. (5), which written out now reads where the very last factor is the x b derivative of F evaluated at x b = 0. The approximation of Ref.
[18] amounts to ignoring the fact that the derivatives of F s are evaluated at a different point in the last two lines, which gives This has the form of a Liouville equation as noticed in Ref.
[18] and so its solution can be written as where the components of the vector functions ψ s and ψ b evolve with τ according to from the initial conditions The function F s at τ = 0, which as before we write without a time argument, is given by the analog of Eq. (12): The corresponding memory as defined in Eq. (19) is FIG. 13. Channel decomposition of the memory terms for the ZMs trajectory shown in Fig. 6. (a) and (b) show memory terms affecting Nkx2.2 and Olig2, respectively. Each color represents a memory channel as indicated. The memory originates from a particular subnetwork species "sending" memory through a specific bulk species; the effect then propagates within the bulk and returns via a specific bulk species (see legend). The two most salient memory functions are Nkx2.2 to Pax6 and then returning through Pax6 into Nkx2.2 (thick yellow line), and Olig2 to Irx3 and then returning through Pax6 into Nkx2.2 (thick blue line). (b) shows a large memory contribution acting on Olig2 via the channel through Irx3 (yellow line). However, in this case, the drift for Olig2 (c) is so large that the relative effect of this memory channel remains nonetheless small. where all R s , R b and the derivatives are evaluated at (x s , 0). Gouasmi et al. propose to find these derivatives numerically, but in fact a closed form expression can be obtained as follows. Applying the chain rule gives Now note that in the final evaluation we always use x b = 0, which from the differential Eqs. (H6) and (H7) implies ψ s = x s , ψ b = 0 for all τ . Hence, in particular ψ b is independent of x s and so ∂ψ b /∂x s = 0. We also have F s (x s , 0) = 0 from Eq. (H9), which implies ∂F s /∂ψ s = 0. Only the last term from Eq. (H11) thus survives, and it remains to find ∂ψ b /∂x b . By differentiating Eq. (H7) for ∂ψ b /∂τ w.r.t. x b , one finds ∂ ∂τ On the r.h.s., a similar term from the variation of ψ s vanishes because it would be proportional to This difference is zero in the final evaluation at x b = 0 (which implies ψ b = 0). For the same reason, the derivatives ∂R b /∂ψ b = ∂R b /∂x b are evaluated at (x s , 0) and constant in time τ . Collecting these derivatives into a matrix k with elements k b b and using that ∂ψ b /∂x b = δ bb (= 1 for b = b and = 0 otherwise) at τ = 0 gives then as the explicit solution of Eq. (H13), and inserting into Eq. (H12) yields where we have used that ∂F s /∂ψ b = ∂F s /∂x b = ∂R s /∂x b ; this derivative and the factor R b are evaluated at (x s , 0) in the approximation from Ref.
[18] for the memory function.
We do not show here how the above memory approximation performs in our test systems because the nature of the approach can lead to fixed points disappearing after projection or new fixed points appearing. We observed both of these effects in numerical evaluations for the bistable switch from Ref. [30].