Nonlinear delay differential equations and their application to modeling biological network motifs

Biological regulatory systems, such as cell signaling networks, nervous systems and ecological webs, consist of complex dynamical interactions among many components. Network motif models focus on small sub-networks to provide quantitative insight into overall behavior. However, such models often overlook time delays either inherent to biological processes or associated with multi-step interactions. Here we systematically examine explicit-delay versions of the most common network motifs via delay differential equation (DDE) models, both analytically and numerically. We find many broadly applicable results, including parameter reduction versus canonical ordinary differential equation (ODE) models, analytical relations for converting between ODE and DDE models, criteria for when delays may be ignored, a complete phase space for autoregulation, universal behaviors of feedforward loops, a unified Hill-function logic framework, and conditions for oscillations and chaos. We conclude that explicit-delay modeling simplifies the phenomenology of many biological networks and may aid in discovering new functional motifs.

is an ODE, and is a corresponding DDE, where x(t − τ) represents the value of x at a time τ units in the past, making the effect of x on the current rate of change of x delayed by a time τ. The time scales and delays are thus explicit (Fig. 1), better capturing dynamics where intermediate steps are not fast, as compared to ODE models, which may fail to capture real delays, hide their effects by oversimplification, or require additional variables and parameters to predict complex phenomena 17,23,32 . Such delay models have been productive in a variety of biological scenarios, such as processes with many intermediate but unimportant steps 33,34 and agestructured populations 34,35 . Network models with delays have also been explored in various contexts [36][37][38][39] . Intercellular signal transduction, for example, requires multiple time-consuming steps that are not negligible 23,32 . With DDEs, multiple steps ("cascades") within a network can be rigorously simplified into a single step with delay (see below), an approach which has been explored in a variety of biological contexts 18,33,34,[40][41][42] . This makes interpretation of the phenomenology simpler than with ODEs and reduces the number of equations and parameters in the model 18,23 . This idea can be visualized by depicting regulatory networks as directed graphs, with nodes representing biological species, and pointed vs. blunt arrows indicating activation vs. repression, respectively ( Fig. 1). DDE models allow a single arrow to faithfully capture many biochemical steps 18,33,43 , expanding the available dynamics in a model with a reduced set of equations and parameters. DDE regulatory models have in fact been used widely to model a range of biological phenomena such as development 17,44 and hematopoiesis 30,31 . For instance, a 1variable ODE such as Eq. (1) can only produce exponential  Fig. 1 Explicit inclusion of time delays in mathematical models of network motifs can reproduce non-delay models using fewer variables and parameters, but can also lead to more complex behavior. a An example genetic regulatory network including genes X, Y, and Z, regulated by one another either positively (activation, pointed arrow) or negatively (repression, blunt arrow). The model on the left incorporates no explicit delay terms, whereas the model on the right incorporates explicit delay terms τ x and τ z . b Ordinary differential equations (ODEs) and delay differential equations (DDEs) corresponding to (a) with regulation strengths α i , removal rates β i and cooperativities n i . c Numerical simulations of equations in (b) for one set of initial conditions using parameters α x = 1, α y = 1.2, α z = 1.2, β x = β y = β z = 1, n x = n y = n z = 2 for the ODE simulation, and α x = 0.5, α z = 0.3529, β x = β z = 0.5, n x = n z = 2 for the DDE simulations (dashed curves: τ x = 0.8, τ z = 0.1; solid curves: τ x = 3, τ z = 4). Note that delays can cause more complex dynamics (e.g., transient oscillations) compared to models in which effects are instantaneous, and where both models give the same long-term steady state behavior. Left: instantaneous effects modeled using ODEs, Right: delayed effects modeled using DDEs. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-21700-8 growth or decay 45,46 . A corresponding DDE, such as Eq. (2), can also oscillate (with amplitude approaching zero or infinity) and, if nonlinear, lead to stable oscillations, bistability and chaos 13,26,[47][48][49] .
A key challenge in using DDEs is their mathematical complexity relative to ODEs 50,51 . For example, while Eq. (1) is 1dimensional because one initial condition (x at t = 0) determines all future behavior, Eq. (2) is infinite-dimensional, because x must be specified at all times − τ ≤ t ≤ 0 to predict a unique solution 50,51 . Concretely, the solution to Eq. (1) can be written as a single exponential, while Eq. (2) generally must be specified as a sum of an infinite number of exponential terms to satisfy initial conditions.
Despite the challenges, much progress has been made in analytical understanding of DDEs 46,52,53 , and numerical methods exist for simulation 54,55 . We thus see an opportunity to use DDEs to recapitulate dynamics found in ODE solutions of network motif behavior with fewer genes and thus fewer modeling parameters and equations (Fig. 1), a type of "modeling simplicity." In this work, we thoroughly examine the most common network motifs 2,6 with explicit delays and present an approachable, step-by-step view of the mathematical analysis (applying established DDE methods 34,51 ) in order to make such delay equations easy to use for biologists and others. We show the essential role of delays in autoregulation, feedforward loops, feedback loops, multiple feedback, and complex networks, as well as instances where delays may be ignored. A reference summary is provided at the end in Table 2. Finally, we discuss how these results can be applied to understanding fundamental design principles of various natural biological systems.

Results
Our results are divided up into 7 sections corresponding to 7 different regulatory networks of increasing complexity. An overall reference table (Table 2) is included in the discussion.
Motif 0: direct Hill regulation. We first describe how we use DDEs to model simple regulation of a gene Y by a gene X as the most basic motif, and then provide a simple yet unified mathematical framework for both activation and inhibition with delay.
Activation and inhibition can both be modeled with a single unified function. We consider a transcription factor x regulating the production of a protein y (Fig. 2). If x activates y, the production rate of y increases with increasing x, generally saturating at a maximum rate α. Often there is an additional cooperativity parameter n which determines how steeply y increases with x in the vicinity of the half-maximal x input value k. In this framework, n = 0 is constitutive production, n = 1 is a Michaelis-Menten regulation, 1 < n < 5 is a typical biological range of cooperativity, and n → ∞ is a step-function regulation. If x represses y, the same conditions hold except that the production rate of y then decreases with increasing x. A standard quantitative model for this behavior is called the Hill function, and serves as a good approximation for many regulatory phenotypes in biology, including transcription and translation rates, phosphorylation, enzymatic activity, neuronal firing, and so forth 56,57 . In general there may also be a leakage rate α 0 that yields constant y production in the absence of x. The concentration of y is also generally removed at a rate β proportional to its concentration. This removal term can represent many biophysical processes, such as degradation, dilution, compartmentalization, or sequestration 58,59 ; for simplicity we mainly use the term "degradation." Together these can be written as: If x instead "represses" y instead of activating it, the form is similar: For biologically meaningful results, all variables and parameters in these equations should be real and non-negative. Note, however, that the activator case is equivalent to the repressor case with n < 0. In this paper we will therefore allow n to be negative, which is simply a notational modification used in order to combine the two cases. With this notation, the effective cooperativity is |n|. Thus, we have: which provides a unified, single-function description for both activators and inhibitors (and constitutive expression), providing a powerful mechanism to analyze both cases simultaneously as we will do throughout this work. The negative value of n allowed here should not be confused with the term "negative cooperativity" used in the binding kinetics   (7)) of activated (left) and repressed (right) Y following a rise in X that exponentially approaches a new steady state (η X = 6) from zero (governed by _ XðTÞ ¼ η X À XðTÞ). Note that the finite value of η X leads to an effective leakage slightly greater than ϵ for the repressor case. For the activator (n = −2), we chose η = 5.6528 to match esthetically between the X and Y steady states, such that η X ¼ ϵ þ η=ð1 þ η n X Þ. For the repressor (n = 2), we similarly chose η = 5.5 to match esthetically between the X steady state and Y initial condition, such that η X = η + ϵ. In both cases, γ = 2, ϵ = 0.5, and η X = 6. Initial conditions were X = Y = Z = 0.01 for all T ≤ 0. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-21700-8 ARTICLE NATURE COMMUNICATIONS | (2021) 12:1788 | https://doi.org/10.1038/s41467-021-21700-8 | www.nature.com/naturecommunications literature 60 , which in our notation would refer strictly to a negative value of n with magnitude less than one (i.e., − 1 < n < 0). The Hill function has an inflection point for all |n| > 1, which allows there to be 3 fixed points for n < 0 when feedback is introduced below ( Supplementary Fig. 1).
Delays may be present in regulation, but are not modeled in removal. Converting Eq. (5) from ODE to DDE to account for explicit delay τ in regulation, we replace x(t) with x(t − τ) in the regulation (Hill) term, but not in the removal term, which is directly dependent on y itself, and thus not expected to have any delay. A delayed removal term can also in general lead to negative values, limiting its use in a biological context. This form of regulation is quite general, but for concreteness we will generally refer to quantities like x (or y) as the concentrations of some protein x. Thus we arrive at an explicit-delay DDE model of activating or repressive biological regulation as follows: αk n k n þ x n ðt À τÞ À βyðtÞ À1 < n < 0 ðactivatorÞ Nondimensionalizing yields 4 key parameters for any regulation.
We can non-dimensionalize Eq. (6) by dividing all concentrations by the half-maximal input concentration k and dividing times by the degradation time 1/β, which has the effect of measuring concentrations in units of k and times in units of 1/β. This yields: with dimensionless variables X = x/k, Y = y/k, T = tβ, and dimensionless parameters ϵ = α 0 /kβ, η = α/kβ, γ = τβ, and n. We thus reduce the number of parameters from 6 to 4, and as discussed below, primarily η and γ are important. In this form, η is a normalized "regulation strength" and γ is a normalized delay. These dimensionless parameters can be considered "big" or "small" by comparing their value relative to unity. For example, γ = τβ~1 implies that the delay (τ) and degradation (1/β) times are approximately equal. In subsequent sections we generally explore parameter values between 0.1 and 10. The dynamics of direct activation and repression via Eq. (7) are demonstrated by numerical simulation in Fig. 2c. For simplicity, here and in all other simulations, we use constant initial conditions (variables constant in time for T ≤ 0) unless otherwise specified.
Motif I: cascade (sequential regulation). A common network motif in many biological networks is the cascade (Fig. 3), a series of regulatory steps (i.e., x regulates y, which regulates z, etc.) 2,6,61 .
Since each regulator must reach the corresponding half-maximal input value k before significantly affecting the next item in the cascade, each step adds an effective delay. The following analysis Multi-delay DDE model Fig. 3 A delayed model of cascades recapitulates the implicit delays in ODE models. a A cascade is a linear sequence of regulation steps, here X regulating Y regulating Z. b Standard models of cascades use ODEs (top), in which each step leads to a characteristic delay in the products Y and Z, which increases with each step based on the half-maximal inputs and degradation rates for each step. In an equivalent DDE model (middle), similar behavior is accomplished by replacing the explicit cascade of implicit delays with a single-step regulation including an explicit delay. A cascade in which each step contains an explicit delay (bottom) behaves analogously to delayed direct regulation (as in middle), with the final step delayed by the sum of delays in each step. η X = 1.5, η Y = η Z = 2.17, β X = β Y = 1, β Z = 0.667, n X = n Y = − 2. For the bottom graph in (b), each step is governed by Eq. (9) with identical parameters. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-21700-8 falls within a range of methods used to replace intermediate steps with a delay 18,41,42 .
Delayed direct regulation approximates cascades of non-delayed regulation. For a non-delayed cascade motif in which X regulates Y, which in turn regulates Z (Fig. 3), one can write down a nondimensional set of governing equations as follows: in which β x , β y , and β z are the dimensional degradation rates of X, Y, and Z, respectively (Supplementary Note 1). We have included here an explicit input function X for concreteness, but more generally we can consider any input X(T) with characteristic timescale 1/β x . If Y changes quickly compared to X (i.e., β x /β y ≪ 1), a first approximation to arriving at a simplified model that ignores Y is to replace the appearance in _ ZðTÞ of Y(T) with the X-dependent pseudo-steady state Y pss (X(T)) found by setting _ YðTÞ ¼ 0. A more precise estimate with the looser demand that ðβ x =β y Þ 2 ( 1 can be found by using the same pseudo-steady state, but delayed in time as Y pss (X(T − γ)), with γ = β z /β y (Supplementary Note 1). Intuitively, this delay takes into account the fact that on the slower timescale of X, Y effectively follows the dynamics of X, but the delay accounts for the non-zero time β z /β y that it takes for Y to respond to such changes in X (the β z deriving merely from normalizing time T = tβ z ). By removing Y we lose certain information, such as the ratio of Z or X to Y, but the inclusion of a delay provides for a reasonably quantitative approximation of the output Z's dynamics as a function of the input X.
After plugging the delayed pseudo-steady state of Y into _ ZðTÞ, one can match the values of the composite Hill-within-Hill function at X = 0, X = 1, and X → ∞ as well its the slope at X = 1 compared to a single Hill function with leakage, and thereby approximate the cascade as a single-step regulation of Z by X (Supplementary Note 1): where the combined regulation parameters are (Supplementary Note 1) The k can be removed by renormalization of X. Equation (10) has several biological consequences. First, the overall Hill coefficient h is negatively proportional to the product of individual cooperativities n X , n Y . This makes the cascade activating for either two activators or two repressors, and repressive otherwise. It also makes the cascade more cooperative than its components, in line with past analyses 62 and as found experimentally 61 . Third, the leakage is zero for n Y < 0 and negligible (i.e., α 0 ≪ α) for n Y > 0 if η Y ≫ 1. This means Y must be produced in sufficient amount during the delay to repress Z. We generally ignore leakage in the main text; for a more detailed discussion see Supplementary Note 3 and Supplementary Fig. 2.
The total delay in a cascade equals the sum of individual delays. For cascades with delays at each step, the same analysis implies that total delay equals the sum of individual delays (β z /β y + γ y + γ z , with γ y the intrinsic delay between X and Y and γ z the intrinsic delay between Y and Z). Figure 3 shows simulations of ODE, equivalent DDE, and multi-delay DDE cascade models, matching the above analytical results.
Motif II: autoregulation. Autoregulation, one of the most common biological motifs 2,63 , describes a single biological species that regulates its own production (Figs. 4a-c).
The complete phase space for autoregulation demonstrates the quantitative and qualitative importance of delays. Based on Eq. (5), the governing equation for such a system with delayed regulation is given by setting Y = X (output equals input) in Eq. (7): This equation has four parameters (η, γ, ϵ, and n). Since leakage must be small relative to regulation (i.e., ϵ/η ≪ 1) for regulation to be strongly effective ("activated" rate much greater than "non-activated" rate), we will focus here on the case with no leakage (ϵ = 0). We treat non-negligible leakage in the supplements (Supplementary Note 3). Note that in general Eq. (11) has no closed-form solution.
Fixed points do not depend on delays. Fixed points X(T) = X * do not change with time, implying _ XðTÞ ¼ 0 and X(T) = X(T − γ) = X * . In Eq. (11), this implies that production must balance degradation (Fig. 4d). For repressors (n > 0) there is only a single fixed point for biologically relevant parameter values. For activators (n < 0), there can be 1, 2, or 3 fixed points (2 fixed points is a border case for −n > 1). Explicitly, the fixed point values are given by or X * = 0 (for activators only). Again, we assumed ϵ = 0 for simplicity. Note that these fixed points only depend on 1/η = βk/α and ϵ/η = α 0 /α, and have no dependence on the delay time τ.
Linearization is sufficient to determine bifurcations in qualitative behavior. Small disturbances from fixed points δX(T) = X(T) − X * ≪ 1 linearize Eq. (11) (Supplementary Note 2). Assuming a solution δXðTÞ ¼ A expðλTÞ yields a transcendental "characteristic equation" for the eigenvalues λ of Eq. (11): The constant function of the fixed point value M(X * ) is defined as whose sign is given by n (positive for repressors and negative for activators).
To determine stability of the fixed points, we solve for conditions in which λ crosses the imaginary axis (Re λ ¼ 0), occuring at either λ = 0 (a saddle-node bifurcation) or λ = iω (a Hopf bifurcation).
We can see that Eq. (15) separates monostability from bistability by noting that the production and removal curves (Fig. 4d) are tangent when Eq. (15) holds (Supplementary Note 2), yielding two stable fixed points when η is above this boundary and only , and it ceases to be a fixed point for n > 0. Note that this boundary does not depend on delay; however, that does not preclude other features such as basins of attraction from depending on delay 64 . Viewed as a function of n, Eq. (15) limits bistability to −n > 1 and the boundary to 1 ≤ η ≤ 2 ( Supplementary Fig. 1b). Since this function is non-monotonic, it is also possible (although biologically perhaps unrealistic except on an evolutionary timescale), to hold a value of 1 < η < 2 and decrease n to a point where the bistability is lost, and then decrease n still further until bistability is gained again.
This boundary has both horizontal (η) and vertical (γ) asymptotes, given as: Thus oscillations require both a minimal regulation strength and a minimum delay. The vertical (γ) asymptote approaches zero as n → ∞, so large regulatory strength can achieve oscillations even with minuscule delay if cooperativity is extremely steep.
The oscillation period can be approximated by linearization near consecutive maxima and minima of the oscillation (Supplementary Note 2). This method yields a period of approximately 2(γ + 1), or 2(τ + 1/β) in dimensionful terms (compare also to 17 ). Biologically, this means that the concentration X is pushed from high to low (or vice versa) after the delay Fig. 4 The complete phase diagram for the autoregulation network motif has analytically derivable parameter regions for bistability, monostablility, monostablility with damped oscillations, and oscillations. a The autoregulation network motif, with a dotted arrow indicating either (b) self-activation or (c) self-repression. The two cases are given by n < 0 (activation) and n > 0 (repression). d Stable (black circles) and unstable (white circles) fixed points for activator and repressor cases are given by the intersection of the regulation (solid) and degradation+leakage (dashed) lines. Note that activators can have 1, 2, or 3 fixed points, whereas repressors always just have one. e Parameter space showing all possible behaviors for autoregulation with delay. Shading shows simulation results (with an interval of 0.1 for both γ and η axes) and curves show the analytically derived bifurcation boundaries. See Supplementary Fig. 3a for cases −3 ≤ n ≤ 3 and Supplementary Fig. 3b for boundaries in η vs. n parameter space. f Representative simulation curves for the four qualitatively different behaviors, with different colors representing different initial conditions (see "Methods"). ϵ = 0. time (τ) has elapsed and the concentration has equilibrated to its new value via degradation (1/β).
Damped oscillations are expected when the largest eigenvalue satisfying Eq. (13) has non-zero imaginary part. This is never true for n < 0, and is true for n > 0 above the curve denoted "spiral boundary" in Fig. 1e and Supplementary Fig. 3. This curve is given as follows: which is derived in Supplementary Note 4. It approaches η = 0 for large delay and has a vertical asymptote at nγe γ+1 = 1, so there is some non-zero delay that has non-oscillating decay for any strength η and finite n.
Putting together the derived boundary curves yields an η-γ phase space for each n, displayed in Fig. 4e and Supplementary Motif III: logic. A general class of functions used to describe natural 2,58,65-67 and synthetic 68,69 biological networks are logic gates, which have two inputs regulating a single output (Fig. 5). Gates exhibit either high ("on") or low ("off") output depending on whether inputs are on or off. For example, the AND gate specifies high output only if both inputs are on. In this section we provide a specific DDE-based framework that covers 14 out of 16 possible 2-input logic operations, and show that these operations form a continuous 2D parameter space.
A two-parameter summing function reproduces all 2-input monotonic logic gates. We first write out nondimensionalized equations corresponding to the logic gate motif as depicted in Fig. 5a. We assume that the degradation constants (β) for Z and R are equal for simplicity, and that there is no leakage.
We describe the regulation of Z by X and Y using a sum of Hill terms. Logic gate behavior can be captured with other approaches such as a product of Hill terms, or summation within a single Hill term 58,66,67,71 , each representing subtly different biology. We choose the separate Hill term approach as it describes many logic functions simply by tuning regulatory strengths, and can be extended to include multiplicative terms (Supplementary Note 5). Caveats include multiple states for Z, requiring additional binarization via R (Fig. 5c), as well as poor response to ratiometric inputs, discussed in the next section on feedforward motifs.
Using Eq. (19), we can characterize the motif logic based on the idea of dynamic range matching 72 . Every regulator in Eq. (19) is effectively compared against unity in the denominator of the Hill function for its corresponding output. For instance, Z provides an "on" or "off" signal to R if Z > 1 or Z < 1 respectively. Z can take on values below 1 as long as η Z1 ( maxðX jn 1 j Þ and η Z2 ( maxðY jn 2 j Þ, otherwise Z will always activate R, as indicated by the areas marked TRUE in Fig. 5. Let us say that X and Y settle on steady-state values X * ≡ η X , Y * ≡ η Y , as inputs to our logic gate. A value of η X or η Y significantly greater than 1 is then "high" (true, 1), and "low" (false, 0) if much less than 1. We then want to determine whether Z * is greater than (true) or less than (false) 1. For example, if n 1,2 > 0 (two repressors), Table 1 gives the possible steady states of Z. If η Z1 and η Z2 are greater than 1, these steady states approximate a NAND gate. If they are less than 1, but sum to greater than 1, the steady states instead approximate a NOR gate.
Equivalent calculations for all variations of n 1 , n 2 , η X , η Y , η Z1 , and η Z2 (Supplementary Table 2) show that logic is dependent on only two parameters, sgn ðn 1 Þη Z1 and sgn ðn 2 Þη Z2 . These form a 2D parameter space that summarize the logic into a convenient chart (Fig. 5b). This implies evolution or other tuning can smoothly convert between logic gates (see also 66,67 ). The chart is symmetric on interchange of the two inputs (η 1 ↔ η 2 , n 1 ↔ n 2 ), and anti-symmetric (that is, application of logical NOT) by conversion of an activator to a repressor or vice versa (n i → −n i ). Fourteen of the 16 possible two-input logic functions are represented. All fourteen are monotonic in that Z(X, The two nonmonotonic gates, XOR and XNOR, are not represented by this simple logic motif, because summation (addition of Hill terms) is monotonic. However, they can be constructed by connecting several of these gates 70,72 .
The logic parameter space can be divided into AND-type, OR-type, and single-input functions. Only the 8 regulatory functions on the diagonals of Fig. 5 (excluding FALSE and TRUE) make use of both inputs. These eight can be further divided into positive or negative regulation on each arm (4 possibilities) in conjunction with an AND or OR gate (8 possibilities total). Specifically, "ORtype" logic applies when both η Z1 > 1 and η Z2 > 1 (OR, NAND, and both IMPLY gates), and "AND-type" logic applies when both η Z1 < 1 and η Z2 < 1 (AND, NOR, and both NIMPLY gates). This can be seen mathematically by using Boolean logic simplification. For example, It is important to note that this logic scheme is an approximation, and in reality the sum of two Hill terms provides a form of fuzzy logic 73 . That is, if the inputs X or Y are close to 1, or if the regulatory strengths η Z1,2 are close to 1, then Z will also be close to 1 for some input combinations (Supplementary Note 5).
Note that if X(T) and Y(T) are independent, they can be timeshifted in Eq. (19) by γ 1 and γ 2 , respectively, showing that the dynamics do not depend on delays. This is not true if X and Y are dependent (e.g., X = Y), which leads to interesting dynamics that we examine below in feedforward loops and double feedback.
Motif IV: feedforward loop. Equipped with a multivariable generalization of 1-variable Hill regulation (Eq. (19)), we turn our attention to the feedforward motif (Fig. 6), a non-cyclic regulation network in which an input, X, regulates an output, Z, via two separate regulation arms. One arm is "direct", in which X regulates Z in a single step, while the second arm is "indirect", with X regulating an intermediate Y, which in turn regulates Z. In this way, the first arm "feeds forward" past the cascade (Fig. 6a-b). The motif is found commonly in biological networks 6,74 , comprising about 30% of three-gene regulatory interactions in transcriptional circuits 75 . Feedforward loops are conventionally described as "incoherent" if X activates Z through one arm but represses through the other, and "coherent" otherwise (Supplementary Table 3). In this section, we show that the essential behaviors of feedforward loops are due to a difference in delays between the two inputs to Z and the logic function between the two inputs to the output Z.
Feedforward dynamics depend only on the difference in time delay between two regulatory arms and their 2-input logic. With a DDE model, we can drop the cascaded intermediate Y, focusing solely on Z as regulated by X via two arms that have unequal delays (Fig. 6a, b). Looking at Fig. 6a, we can then write down a simple equation for the regulation of Z, by setting Y = X in Eq. (19) (Supplementary Note 6): Because each regulation term is controlled by the same input (X), we cannot in general normalize the half-maximal input to 1 in both terms as we did for logic (Eq. (19)); K = k 1 /k 2 represents the ratio of the two input scales. Likewise, we cannot entirely time-shift out the delays; instead, shifting X backwards in time by γ 1 (XðTÞ ¼ XðT þ γ 1 Þ), we see that the behavior of the feedforward equation actually depends only on the difference in regulatory delays, Δγ = γ 2 − γ 1 rather than each delay independently: For clarity we will from now on refer toX as X without the hat unless specified.
To understand the behavior of Eq. (21), we will analyze the output of Z to two types of input X: a step input, and a continuous oscillation. Termed the step and frequency responses, respectively, these are often used in control theory to characterize input-output relations 76 .
Step-pulse response of feedforward motifs can be solved analytically and demonstrates filtering and pulse-generation behaviors. We considered a step input (Fig. 6c) to include both a step on (from low to high X) as well as a step off (from high to low X): where the Heaviside function is defined such that Θ(T) = 0 for T < 0 and Θ(T) = 1 for T > 0. Eq. (22) represents a square input pulse of width ω starting at X 0 and reaching to η X . This is an onpulse if η X > X 0 and an off-pulse if η X < X 0 .
The lack of feedback in feedforward loops and the fact that the square-pulse input of Eq. f ðTÞ ¼ 1 À e ÀT À Á ΘðTÞ À 1 À e ÀðTÀωÞ ΘðT À ωÞ Here, Z ss is the steady-state value of Z, Z n κ are the magnitudes of deviation from steady state due to each arm, and f(T) specifies where the responses from the two regulation arms are active. In terms of the original X (as opposed to the time-shiftedX), the output in Eq. (24) is shifted to the right in time by γ 1 relative to the original input X(T).
There are several interesting results to note in this equation. First, the cooperativities (n 1 , n 2 ) and thresholds (K) affect the response magnitude, but not the dynamics (f), because the input signal is already infinitely steep.
Second, the η values only affect the steady-state Z ss and the response logic (as sgn ðn i Þη i , see logic discussion). Only 8 logic gates correspond to feedforward loops, in that they respond to both inputs: the four AND-type functions and four OR-type functions.
Third, f indicates where responses are active, varying from 0 far outside a pulse, aprroaching 1 inside for wide pulses (ω ≫ 0). The output then has four responses (pulse start and pulse end for each arm), in which Z moves away from Z ss by as much as Z n κ (pulse on) before returning (pulse off). Responses must be driven externally, unlike autoregulation.
For coherent feedforward motifs, input pulses must be sufficiently wide (ω ≫ 1) to fully activate the output R. This implies that coherent feedforward motifs act as filters against short signals, either short on signals (AND and NAND) or short off signals (OR and NOR). The delay difference Δγ determines how quickly the response starts (AND and NOR logic) or returns to baseline once the pulse is over (for OR and NAND logic), so that large Δγ maintains a response long after the input pulse has subsided.
For incoherent feedforward motifs, the output R is only fully activated if the pulse is wide (ω ≫ 1) and the delay difference is also wide (Δγ ≫ 1), in which case a pulse of output is formed. The output pulse is formed on either the on-step (slow arm IMPLY fast arm, and slow arm NIMPLY fast arm) or the off-step (fast arm IMPLY slow arm, and fast arm NIMPLY slow arm). This means that incoherent feedforward motifs act as pulse generators. The delay difference Δγ determines the minimum width of input pulse that will produce a response in R.
Many of the feedforward behaviors described in this section corroborate earlier results from ODE models 2,2,58,74,75,77,77,78 , but are derived here with just a single equation. In particular, it is well known that coherent feedforward loops have filter capabilities and that Fig. 5 A simple approximation for digital logic using a sum of Hill terms recapitulates all monotonic logic functions in a single parameter space. a A prototypical regulatory network involving logic where X and Y both regulate Z, which must integrate the two signals using some logic before it can in turn activate a downstream reporter R. b Parameter space showing regions where regulation approximately follows 14 of the 16 possible 2-input logic functions depending on the strength of two single-variable Hill regulation terms (η Z1 : regulation of Z by X, η Z2 : regulation of Z by Y). Network logic can be smoothly altered by varying the parameters (η Z1 , η Z2 ), with a change of sign in (n 1 , n 2 ) required to switch quadrants. The bottom-left quadrant shows that very weak regulation in both terms leads to an always-off (FALSE) function, weak regulation in one arm only leads to single-input (X, Y) functions, strong regulation in both arms leads to an OR function, and regulation too weak in either arm alone to activate an output but strong enough in sum leads to an AND function. The other three quadrants are related by applying NOT to one or both inputs, with function names related by de Morgan's law 70 NOT(X OR Y) = NOT X AND NOT Y. In particular, X IMPLY Y = NOT(X) OR Y, X NIMPLY Y = X AND NOT(Y), X NOR Y = NOT X AND NOT Y, and X NAND Y = NOT X OR NOT Y. Truth tables for all 16 logic gates are provided in Supplementary Table 1 for reference. The two non-monotonic logic functions, X XOR Y and X XNOR Y, are those 2 of 16 not reproduced directly using this summing approximation. They can be produced by layering, e.g., NAND gates 70 . c Representative time traces for AND (η Z1 = η Z2 = 0.9) and OR (η Z1 = η Z2 = 1.8) gates with n 1 = n 2 = −2, n 3 = −20, η R = η Z1 + η Z2 . The function sgn ðnÞ ¼ þ1 when n > 0, sgn ðnÞ ¼ À1 when n < 0.  Frequency response of feedforward motifs can be solved analytically and demonstrates low-and band-pass filtering capabilities. Biological regulatory networks often encode information as the change in frequency of an oscillating input 80 , which has been suggested to be more robust to noise than encoding information in absolute concentration [81][82][83] . That feedforward loops filter short square pulses suggests more general frequency-filtering capabilities. We therefore analyze feedforward frequency response to sinusoidal input ( Fig. 6d-g). We show that the response follows the outline of a universal transfer function curve, independent of the logic or delay difference. Instead of the step input analyzed above, here we consider a sinusoidal input which oscillates between zero and 2A (twice the amplitude) at a frequency f > 0. Taking the Fourier decomposition of each Hillregulated term 1/[1 + X n (T)] and plugging into the governing Eq. (21) (Supplementary Note 6) provides the output Z(T) in terms of its magnitude I k and phase ϕ k as a function of frequency: where a x;n k are Fourier coefficients of Hill-regulated terms 1/[1 + X n (T)] at integer multiples k of the fundamental frequency f: for all f > 0. Note that the magnitudes are symmetric to interchange of the two regulation arms (η 1 ↔ η 2 , n 1 ↔ n 2 ), while the phases are not. Looking at frequencies that are integer multiples of 1/Δγ, we can see that both the magnitudes and phases follow a universal envelope not dependent on delays or logic (i.e., no dependence on γ or η values): where (from Eq. (26), at k = 0 or f → 0) I 0 ¼ η 1 a x;n 1 k þ η 2 a x;n 2 k and ϕ 0 ¼ πðsgn I 0 þ 1Þ=2. The actual magnitudes and phases are then modulated from these envelopes with period 1/kΔγ. The fundamental-frequency (k = 1) response (the largest magnitude, see Fig. 6d, e) are shown compared to the envelopes in Fig. 6f (known as Bode plots 76 ).
At frequencies f that are integer multiples of 1/kΔγ, the cosine in I k equals ±1. Because a x;Àn k ¼ Àa x;n k (see Eq. (27), Supplementary Note 6), the output magnitude I k then increases to a local maximum for coherent feedforward motifs (sgn ðn 1 Þ ¼ sgn ðn 2 Þ) and decreases to a local minimum for incoherent feedforward motifs (sgn ðn 1 Þ ¼ Àsgn ðn 2 Þ). The opposite holds for frequencies that are half-integer multiples of 1/kΔγ. For the special case of perfectly balanced incoherent feedforward motifs (η 1 = η 2 , n 1 = −n 2 , K = 1), the magnitudes I k decrease to zero for frequencies that are half-integer multiple of 1/ Δγ; otherwise, the maxima (for coherent) and minima (for incoherent) equal I k (f = 0).
Incredibly, the envelopes (relative to f = 0) have zero parameters, making them universal to all feedforward loops. The frequency at which the magnitude decreases to one-half its value at f = 0 occurs at f ¼ ffiffi ffi 3 p =2πk, with a phase shift of −π/3. Despite the non-dependence on Δγ, the absolute maximum magnitude obtained for incoherent feedforward motifs does depend on Δγ, as does the frequency at which this maximum response occurs. The dependence on delay difference is particularly strong near Δγ = 1 (Fig. 6g). This is because for Δγ ≫ 1, the envelope decreases by half when f is many multiples of Δγ, and thus past many local maxima. The local maximum then occurs at 0 < f ≪ 1 with a maximum above I k (f = 0). For Δγ ≪ 1, on the other hand, the envelope decreases by half before a single sinusoid is complete, so the maximum is at f → 0.
At Δγ ≈ 1, the first sinusoid reaches its maximum while the envelope is descending most steeply, and the maximum occurs at f ≈ 1 with a value highly dependent on Δγ. To activate R, X must go above 1 for at least part of the output oscillation, meaning that I 1 must go above 1 − 〈I〉 = 1 − I 0 . Looking at Fig. 6g, we see that for large Δγ ≫ 1, this occurs for small frequencies, and the feedforward loop acts as a low-pass filter. However, there is a window of Δγ ≈ 1 in which 1 − I 0 is large, at modest frequencies (0.1 < f < 10), corresponding to the hump in magnitude, e.g., around f = 0.5 for Δγ = 1 in Fig. 6f. The incoherent feedforward motif thus acts as a bandpass filter for Δγ ≈ 1, with ineffective regulation outside of modest frequencies.
For coherent feedforward motifs, the frequency at which the first minimum in output response occurs follows the same logic, but the maximum response is always found at f → 0. Since the maximum response is at low frequencies, coherent feedforward motifs act as low-pass filters, as has been shown to be a basic expected behavior of chemical networks 84 .
Motif V: multi-component feedback. The common twocomponent feedback motif 58 is the equivalent of a cascade with feedback (Fig. 7a). We can write down the governing equations for such a system, using the same normalizations as before (see Supplementary Note 7), assuming for simplicity zero leakage and that x and y have the same degradation rates β. This is equivalent to the 3-step cascade Eq. (8) with equal degradation rates and Fig. 6 The feedforward network motif owes its primary functions to a difference in regulatory delays. a The feedforward motif with delays, in which a single output Z is controlled by an input X via two regulatory arms with differing delays. The straight, short arrow represents the "direct arm" with delay γ 1 and the longer, curved arrow represents the "indirect arm" with delay γ 2 > γ 1 . b The ODE model for an incoherent (type 1) feedforward motif, one of 8 possible networks in which the intermediate gene Y is modeled explicitly in the indirect arm. c Simulations of the four feedforward motifs with AND-type and OR-type logic (Fig. 5 and Supplementary Table 1) in response to short and long gain and loss of input signal. Blue curves: inputs (X), orange curves: reporter R activated with high cooperativity by Z. Note that the bottom two rows demonstrate pulse generation, while the top two rows filter short signals. d Response of an incoherent feedforward motif to oscillatory input after initial transients have died away. Z 3 is a 3-frequency Fourier approximation of Z (see E). e Fourier decomposition of Z from (d) by Eq. (26) and by a numerical fit to the data in (d). f Frequency scan (Bode plot) of (d) for 3 values of Δγ, with the theoretical envelopes from Eq. (28). g The maximum amplitude of the motif in (d) over a range of Δγ and the corresponding frequencies at which the maxima occur. Z goes above 1 (activation threshold for R) for a small range of Δγ. For (c-f), η 1 = 0.9, η 2 = 0.7, n 1 = 2, n 2 = −2, n 3 = −20, η R = 2, A = 1.
setting Z = X to form the loop.
Saddle-node bifurcations determine bistability modes. For saddlenode bifurcations, at which bistability begin, we set λ = 0. This can only be satisfied by the Λ + version of Eq. (31) if n 1 < 0, n 2 < 0 (two activators), and by the Λ − version if n 1 > 0, n 2 > 0 (two repressors). A feedback loop of an odd number of repressors cannot have a saddle-node bifurcation. As for linear cascades (Fig. 3), the overall autoregulation is repressive only if it contains an odd number of repressors. The bistability boundary for both cases is then given by (Supplementary Note 7): where f = 1 provides an implicit monotonic mapping between X and Y that restricts η 1 (X) with respect to η 2 (Y). The phase difference ϕ = 0, implying that two-activator loops show bistability with X and Y either both high or both low, because the components of v + are of like sign, whereas tworepressor loops (v − ) show bistability with X high and Y low or vice versa (Fig. 7b, d).
Hopf bifurcations determine transient oscillations in modes restricted by bistability. For Hopf bifurcations, at which oscillations begin, we set λ = iω. The boundaries are given for both eigenmodes by (Supplementary Note 7) where again f = 1 provides an implicit monotonic mapping η 1 (X) and η 2 (Y). These curves lie above the bistability boundary, so oscillations in v + are not observed for two activators nor in v − for two repressors.
Based on Eq. (33), we can see that the primary Hopf bifurcations are equivalent for differing values of the delays as long as the average delay is equal, because the oscillation frequency ω of the eigenmodes only depends (implicitly) on the average delay: In particular, oscillations are possible for activator-repressor loops even when one delay is zero, as long as the sum of delays is greater than zero. On the other hand, the phase difference between the oscillations depends only on the difference in delays given the value ω: Together with the additional phase difference of π for v − , this implies that for equal delays, there are synchronous oscillations for two-repressor loops, anti-synchronous oscillations for twoactivator loops, and π/2-shifted oscillations for activator-repressor loops. This phase relation also holds off the bifurcation boundary, where λ = μ + iω, as ϕ does not depend on μ (Supplementary Note 7). Note also that the overall phase in the expressions for γ are 0 for two repressors, π for two activators, and π/2 for activator-repressor loops. Finally, one can also conclude from Eq. (30) that (Supplementary Note 7) which holds for both the real and imaginary parts of λ. This implies that the oscillations grow most quickly and have the highest frequency when the delays are nearly equal (γ 1 = γ 2 ). This consequently makes the transients most noticeable for equal delays. All these results can be seen in the simulated data (Fig. 7b,d).
Motif VI: double feedback. Multiple feedback (Fig. 8) has been described to lead to interesting dynamics, such as excitablity 85 , quasiperiodicity 71 , and chaos 71 . A full analysis of this motif is beyond the scope of this work, for reasons described below. We focus on chaos, a dynamical behavior characterized by irregular, unpredictable oscillations 45 . It has been suggested that biological systems avoid chaos 86,87 , but also that chaos can describe pathological dynamics such as irregular breathing 26 and epilepsy 27 . There may also be cases where chaos is important to biological function [88][89][90][91] . Mathematical models with 71,92,93 and without 86,94 delays can be chaotic.
Chaotic behaviors are prevented in monotonic regulation with linear or cyclic network topology. Due to cyclic network topology (no more than one loop) and monotonic regulation, none of the motifs so far can yield chaotic dynamics 46 . Either non-monotonic feedback (as in the Mackey-Glass equation 26 ) or non-cyclic topology 71 is required for chaos. Double feedback (Fig. 8a) is a minimal motif fulfilling both requirements, although the nonmonotonicity is sufficient (Supplementary Note 8).
The governing equation for such a double feedback motif is thus given by setting the output and both inputs of the logic equation (Eq. (19)) to X (and letting K = 1 for simplicity): Chaotic behavior is possible for non-monotonic feedback with multiple feedback and disparate delays. Simulation of Eq. (37) demonstrates chaos for some parameters of positive/negative mixed feedback, as evidenced by sustained oscillations with variable maxima (Fig. 8b), occupying an apparently fractal region in phase space (Fig. 8c), and a large number of peaks in frequency space (Fig. 8d). To confirm that the dynamics are in fact chaotic, we calculated the dominant Lyapunov exponent using the Wolf method 95 (see "Methods"), finding a positive (i.e., chaotic) value of 0.0040 ± 0.00055 bits (mean ± standard deviation). We also calculated the box dimension of a reconstructed phase space with coordinates (X(T), X(T − 10), X(T − 20)) used by the Wolf method, yielding fractional dimension~1.81, with 95% confidence interval (1.76, 1.87). These results indicate chaos (see Supplementary Fig. 4).
Varying γ 2 relative to γ 1 = 1 (Fig. 8e) shows regions of presumably chaotic dynamics interspersed with simple oscillation, as found in many chaotic systems 45 . In particular, we only find chaos if the delays in the two arms differ substantially. Both delays appear to be important, not only their ratio or difference, since no oscillations exist when both delays are less than one. A full exploration of the two-delay parameter space is left to future work.
In constrast, we found quasiperiodic dynamics for dual negative feedback ( Supplementary Fig. 5). Qualitatively, the time dynamics appear less pulsatile ( Supplementary Fig. 5b) than dual positive/negative feedback, in line with reports that dual positive/ negative feedback demonstrates excitability 85 . The attractor also has a somewhat different profile ( Supplementary Fig. 5c) that in reconstructed space (X(T), X(T − 10), X(T − 20)) traces out a torus with box dimension~2 and near-zero Lyapunov exponent (see Supplementary Fig. 4), indicating quasiperiodic rather than chaotic dynamics. The Fourier spectrum also exhibits a less dense set of peaks ( Supplementary Fig. 5d).
Many parameter choices can prevent chaos in dual feedback. Double autoactivation (n 1 < 0, n 2 < 0) prevents chaos via a stable fixed point (X = 0) 45 . Identical regulation strength and cooperativity (η 1 = η 2 , n 1 = n 2 ) behaves like autoregulation close to fixed points, with γ = γ 1 + γ 2 , as do identical delays (γ 1 = γ 2 ), with ηM = η 1 M 1 + η 2 M 2 (Supplementary Note 8). Thus, for chaotic solutions to occur via double feedback, there must be at least one negative feedback arm, the delays must differ, and either the cooperativities or strengths must differ. If biological systems evolved to avoid chaos 86,87 , it may be that these conditions are selected against, even if the double-feedback motif is not.
The complete parameter space for double feedback is likely to be highly complex, with boundaries possibly formed via an infinite number of bifurcation curves 48 . In particular, while we demonstrated chaotic and quasiperiodic dynamics for specific parameter values, a more in-depth exploration is warranted to fully understand the dynamics for all parameters. For further exploration of double negative feedback and more complex delayed network motifs that exhibit chaotic dynamics, see Suzuki et al. 71 .

Discussion
In this work, we systematically studied the most common network motifs ( Table 2) by use of DDEs. Compared to the more widely used mechanistic ODE descriptions, we demonstrated numerous instances where inclusion of explicit delays simplifies complex systems into smaller motifs while continuing to capture their key dynamics that are otherwise lost in simplifying complex systems as small motifs. To further modeling simplicity, we also identified unifying descriptions for activators and repressors (Eq. (7)) and for boolean Hill function logic (Fig. 5) that are useful for both DDE and ODE models.
There are multiple findings to highlight. First, multistep cascades with and without feedback reduce effectively to direct, delayed regulation, providing an intuitive interpretation of biological delays as resulting from multi-step processes, and establishing analytical relationships for converting between parameters of the corresponding ODE and DDE models. Second, delays determine many key motif properties, such as the oscillation period in negative autoregulation, and pulse width and frequency cutoff in feedforward motifs. Interplay of multiple delays plays a similar role in multi-component feedback for determining absolute and relative oscillation periods, and in multiple feedback for determining chaotic dynamics. Third, we showed quantitatively not only when delays are crucial to behavior, as in oscillations, but also when they may be safely ignored. Delays do not affect steady states or logic of independent inputs, only a difference of delays between regulation arms is important in determining feedforward behavior, and only the sum of delays determines the onset of transient oscillations in feedback loops.
Using DDE models to contract a network into a much smaller network with equivalent topology (replacing all cascades with delayed, direct regulation) may aid in discovery of large-scale motifs that have important functions or that might otherwise be perceived as statistically insignificant. The motifs in this paper were originally discovered by scanning biological networks for subnetworks of N nodes that occur more frequently than expected by chance 6,11 (this in fact being the original definition of network motifs). However, this discovery method becomes increasingly difficult for N ≳ 5 due to the combinatorial scaling of the number of possible motifs 96 . Performing a similar search on contracted networks would be equivalent to searching the original network for larger motifs, and places a stronger emphasis on topology than on number of components involved (e.g., X → Z being equivalent to X → Y → Z). While this can be done without DDEs, the introduction of delays allows one to perform the contraction while maintaining key information about the original dynamics of retained components.
While our work has focused on the most basic network motif topologies (Table 2), the same techniques can be applied to other more complex networks. Expanding on the eigenmode analysis we performed for multi-component feedback loops, one can describe the linear behavior of an arbitrary network near its fixed points by a matrix form of the characteristic equation Ja ¼ 0 with J ij ¼ ðλ þ 1Þδ ij þ η j M j e Àγ ij λ (Supplementary note 9, Supplementary Fig. 6). This approach allows simple models to provide detailed predictions of complex biological phenomena such as pattern formation 23 .
To clarify the basic role of delays in network motifs, we ignored several phenomena which we suggest that future work explore. These include: (1) Noise, which is intrinsic to many biological systems 97 and known to affect dynamics of delay equations 15,92,98,99 . (2) Time-varying and stochastic delays 42,47,98,100 . (3) Non-constant initial conditions ("histories") for T < 0; for example, autorepression shows in-phase and anti-phase locking with constant histories, while sine-wave histories with randomized phase do not (Supplementary Fig. 7). (4) Non-Hill regulation and complex degradation functions 26,78 , such as zeroth-order 101 , nonlinear 102 , or delayed removal, as well as diffusion or spatial effects 103 . (5) Feedback with multiple delays, which we only briefly analyzed in the double feedback motif, and likely have very complicated parameter spaces 48,49 .
Overall, we believe that our work may help resolve fundamental biological and engineering questions regarding a variety of phenomena, including transcription factor networks 6,17,77 , cell cycles 87 , other biological clocks 14,18,25 , and pattern formation 23,44 . Multiple feedback analysis may determine whether chaos affects the cell cycle 87 , or whether biology avoids possibly chaotic motifs. Delays in multicellular signaling 23,32 may distinguish among models of pattern formation 16,24,32,103,104 . Finally, exploring "contracted" networks with delays may uncover entirely new functional network motifs that are larger and more complex than currently known 96 .

Methods
Analytics and numerical simulation. Analytics were in general performed by hand, and checked for validity using Mathematica. Numerical simulations were run in Matlab using the dde23 delay differential equation solver for DDEs and ode45 for ODEs. Simulating activators as repressors with n < 0 technically fails when x is identically zero (Eq. (5)), since that would imply division by zero, but the limit as x goes to zero causes the regulation term to be zero, which is the same result as assumed by our notation. An initial value of exactly zero for x can thus lead to a divide-by-zero error in simulations, and so initial conditions of exactly zero were not used, as that case is an uninteresting fixed point for activators in any case. Note also that the consitutive case for Eq. (5) is degenerate, in that n = 0, α ≠ 0 is equivalent to n ≠ 0, α = 0 with α 0 → α 0 + α/2.
Phase plot simulations and analysis. For autoregulation phase plots, simulations were run with 100 constant-history initial conditions spread logarithmically between 10 −4 and 2η and run from T = 0 to T = 100(γ + 1). Solutions were considered stable if for all 100 simulations the maximum absolute value of the discrete derivative in the last three-quarters of the simulation time was less than 0.1. Stable solutions were sub-categorized as bistable if a histogram of final values over all 100 solutions had more than 1 peak. Solutions were considered oscillatory  Figure 3.
II: Autoregulation Eq. (11) Full parameter space derived. Delays matter only for negative autoregulation. Figure 4, Supplementary  Fig. 3, Supplementary Fig. 1, Supplementary Fig. 2 Capable of pulsing, signal filtering by input pulse width or frequency. All signal processing behavior is due to logic and difference in delays between arms. Figure 6.
V: Feedback loop Eq. (29) Delay sum governs presence of oscillations, which are transient for two repressors (synchronous) and two activators (anti-synchronous). Delay difference governs phase. Full parameter space derived. Figure 7.
if the average Fourier transform of the last three-quarters of the simulation time for all 100 solutions had more than zero peaks with amplitude (square root of power) greater than 100. Solutions were considered spiral if this oscillation condition held for the first one-quarter of the simulation time only. For two-component loops, initial conditions were used that ranged between 0 and maxðη 1 ; η 2 Þ, for equal X and Y and for apposing X and Y. Bistability was determined as for autoregulation, and a cutoff of 0.05 was used to determine "low" values. All simulation histories were constant except where indicated in Supplementary Fig. 7. Specific parameter values and simulation details are given in the figures and/or made explicit in the MATLAB code in Supplementary Data 1.
Lyapunov exponents and box dimensions. For calculating dominant Lyapunov exponents, we used the Wolf method 95 with recommended parameters from the MATLAB script "Wolf Lyapunov exponent estimation from a time series" (version 1.2.0.1) provided by the authors on the Mathworks FileExchange (https://www. mathworks.com/matlabcentral/fileexchange/48084-wolf-lyapunov-exponentestimation-from-a-time-series). In particular, we used an embedding dimension of 3 and a phase space reconstruction delay of 10 for a dataset with~350 orbits of 50-60 data points per orbit. The last 50 iterations of the algorithm were used to generate a mean and standard deviation of the estimated dominant Lyapunov exponent. A full set of parameters can be found in the included code. For calculating attractor dimensions, we used the box-counting method, whose code is also provided in Supplementary Data 1, on the reconstructed phase space (X(T), X(T − 10), X(T − 20)) used for calculating Lyapunov exponents. A linear regression was performed using MATLAB's fit function on half the linear domain to generate a mean and confidence interval (using MATLAB's confint) for the slope (i.e., box dimension) between numbers of boxes covering the attractor versus length of each box.