Cumulative dose responses for adapting biological systems

TL;DR

This paper introduces a new way to study biological adaptation by analyzing cumulative dose responses over time, revealing insights into how systems adapt to inputs.

Contribution

The paper introduces the concept of cumulative dose response (cDR) and shows that only negative integral feedback motifs can produce non-monotonic cDRs.

Findings

Incoherent feedforward loop motifs yield monotonic cDRs and are inconsistent with experimental data.

Negative integral feedback motifs can produce non-monotonic cDRs, aligning with observed cytokine accumulation.

Cumulative dose response provides a new framework for analyzing adaptation in biological systems.

Abstract

Physiological adaptation is a fundamental property of biological systems across all levels of organization, ensuring survival and proper function. Adaptation is typically formulated as an asymptotic property of the dose response (DR), defined as the level of a response variable with respect to an input parameter. In pharmacology, the input could be a drug concentration; in immunology, it might correspond to an antigen level. In contrast to the DR, this paper develops the concept of a transient, finite-time, cumulative dose response (cDR), which is obtained by integrating the response variable over a fixed time interval and viewing that integral—area under the curve—as a function of the input parameter. This study is motivated by experimental observations of cytokine accumulation under T-cell stimulation, which exhibit a non-monotonic cDR. It is known from the systems biology literature that only two types of network motifs, incoherent feedforward loops and negative integral feedback (IFB) mechanisms, can generate adaptation. Three paradigmatic such motifs—two types of incoherent loops and one integral feedback—have been the focus of much study. Surprisingly, it is shown here that these two incoherent feedforward loop motifs—despite their capacity for non-monotonic DR—always yield a monotonic cDR, and are therefore inconsistent with these experimental data. On the other hand, this work reveals that the IFB motif is indeed capable of producing a non-monotonic cDR, and is thus consistent with these data.

Figures

Response functions y(t) plotted against time. As an illustration, three outputs y1(t), y2(t) and y3(t) are shown, corresponding to three inputs u1, u2 and u3, respectively, and their values at a specified time T are shown on the vertical axis. In an adaptive system, the values of all the yi(t)'s converge to the same value y^ when t→∞. For each yi(t), solid curves are used for behaviours until the specified time t=T, and dashed lines for the continuation to t=∞. The area of the shaded region represents the integral ∫0Ty1(t)dt.

Left: Dose response (DR) at time T, obtained from time-resolved data as in figure 1. Again, the vertical axis represents the evaluation of responses y(t) at a specific time t=T, but now these values are plotted against the respective input u, instead of against time. In this instance, DR is not a monotonic function; for example, u1<u2 but y1(T)>y2(T). Right: Cumulative dose response (cDR) at time T; now the vertical axis shows the integral (area under the curve) z(t)=∫0Ty(t)dt of the response, again plotted against the input u. For example, z1(t) is the area shaded in blue in figure 1. This particular cDR is monotonic. For example, u1<u2<u3 and z1(T)<z2(T)<z3(t), because in figure 1 the area under the red curve y2(t) is larger than the area under the blue curve y1(t), even though the final value y2(T) is smaller than y1(T).

Cumulative dose response based on average of three experiments. Plot uses experimental data from [6] (see also fig. 2B in that paper). Horizontal axis denotes concentrations of the input (in units of ligand in ng/well).

Plots of DR (y(t)) and cDR (∫0Ty(t)dt) for the example in equations (1.1) and (1.2) . The initial conditions are x(0)=0, y(0)=1/10, and the time horizon is T=1.5. Using logarithmic scale on inputs, for comparison with experimental plots. Observe that the DR is non-monotonic, yet, surprisingly, the cDR is monotonic. Our main theorem proves monotonicity of the cDR in general, for all IFFL1 and IFFL2 systems.

Plot of cDR (∫0Ty(t)dt) for the IFB example in equations (1.3) and (1.4). The initial conditions are x(0)=0.1, y(0)=6, and time horizons shown are T=3,4,5,6. Observe that, just as with the experimental data plotted in figure 3, the cDR is more monotonic (on the shown ranges, at least) for smaller time horizons T. Using logarithmic scale on inputs, for comparison with experimental plots.

Top: cDR plots of individual experiments, and measured at different times. Bottom: adaptation behaviour in individual experiments. Output y(t) is estimated from individual cumulative z(t) plots in respective top panels.

Three examples of systems: (a) is an IFFL with degradation enhancement, (b) is an IFFL with production inhibition and (d) is an IFB system. The input u is assumed to be a positive constant, and x, y are abundances of a quantity of interest, such as concentration of a protein or mRNA. Arrows ‘→’ indicate positive effect (activation) and blunt edges ‘⊣’ denote negative effects (inhibition). Both (a) and (b) have the same qualitative relation, schematically represented by diagram in (c), between activation of x and y by the input u, and inhibition of y by x, but their dynamical properties are very different. Dynamics can be described by pairs of differential equations for the abundances x(t) and y(t) as a function of time, as shown below the diagrams. State variables x(t) and y(t) are taken to be non-negative. When K=0, it is assumed that x(t)>0 for all t. In these equations, α, β, γ, δ and K are all positive constants.

Plot of the map t↦u∫0t∂uy(s)ds for three values of u. Note that for each u, we have u∫0t∂uy(s)ds≥0 for all t≥Tu.

Plot of the maps u↦S(u) and u↦T(u) with x0=p=1.0.

Plot of the map u↦∫0Ty~u(t)dt for some values of p and T. Note that this map exhibits non-monotonicity in these cases.

Equations

Full text

Introduction

The capability to adapt and to formulate appropriate responses to environmental cues is a key factor for the survival of life, at every level from individual cells, to organisms, to societies [1–4]. A delicate balance is needed in this process: organisms maintain tightly regulated levels of vital quantities, even in the face of variations to be counteracted upon, a property sometimes called homeostasis or adaptation, all the while being able to detect and react to changes in the environment. Underlying adaptation at the cellular level are dynamical signal transduction and gene regulatory networks that measure and process external and internal chemical, mechanical and physical conditions (ligand, nutrient, oxygen concentrations, pressure, light, and temperature) eventually leading to changes in metabolism, gene expression, cell division, motility, and other characteristics. These mechanisms enable organisms to display transient responses that gradually return to a baseline activity level when presented with relatively constant input stimulation, a phenomenon usually called ‘perfect’ or ‘exact’ adaptation [5].

Dose response and cumulative dose response

1.1.

In this article, we continue the study of adaptation mechanisms, with an emphasis on monotonicity properties of an output or reporter variable as a function of an input. We are concerned with ordinary differential equation (ODE) systems that model the interactions of several species, and where there is an ‘input’ (which might represent the dose of a drug or of a genetic inducer), whose level is quantified by the variable $[eqn]$ , and an ‘output’ that is time dependent, and whose magnitude at time $[eqn]$ is represented by $[eqn]$ . The input $[eqn]$ will be assumed to be constant, and we write $[eqn]$ to highlight the dependence of the output on both the input and time. Figure 1 shows three typical responses that one might observe experimentally (in the figure, we write $[eqn]$ instead of $[eqn]$ , in order to simplify notations).

Response functions y(t) plotted against time. As an illustration, three outputs y1(t), y2(t) and y3(t) are shown, corresponding to three inputs u1, u2 and u3, respectively, and their values at a specified time T are shown on the vertical axis. In an adaptive system, the values of all the yi(t)'s converge to the same value y^ when t→∞. For each yi(t), solid curves are used for behaviours until the specified time t=T, and dashed lines for the continuation to t=∞. The area of the shaded region represents the integral ∫0Ty1(t)dt.

The dependence on the initial state will not be indicated explicitly; the initial values of all the species will be fixed at values to be discussed.

One defines (‘perfect’) adaptation to constant inputs as the property that, no matter what is the actual value of the input, the numerical value of $[eqn]$ for large $[eqn]$ is the same, that is, $[eqn]$ for some fixed value $[eqn]$ which does not depend on the particular input $[eqn]$ (see an illustration in figure 1). This value represents a habituated or no-response state, as one achieves when presented with a constant background noise or level of light. In engineering terms, an adapting system is a ‘high pass filter’ that essentially acts on a derivative of the input. Adaptation, by definition, is an asymptotic property, since it ignores finite-time behaviour. On the other hand, transient behaviours, particularly how $[eqn]$ varies with $[eqn]$ at a fixed time $[eqn]$ , are often of interest (What is bethe concentration of a drug in a tumour microenvironment, after 1 h, as a function of the drug dose? What is the size of a tumour after 60 days of the start of therapy, as function of the drug dose? What is the size of the pool of infected individuals in an epidemic population model, as a function of transmission parameters?).

We will call $[eqn]$ , viewed as a function of $[eqn]$ the dose response (DR) of the given system, and denote it as $[eqn]$ . One may perform experiments, exciting a system with different input values $[eqn]$ and measure $[eqn]$ as the final output value, thus obtaining a plot of $[eqn]$ versus $[eqn]$ . The left panel of figure 2 shows several DRs obtained from time-resolved data such as presented in figure 1.

Left: Dose response (DR) at time T, obtained from time-resolved data as in figure 1. Again, the vertical axis represents the evaluation of responses y(t) at a specific time t=T, but now these values are plotted against the respective input u, instead of against time. In this instance, DR is not a monotonic function; for example, u1<u2 but y1(T)>y2(T). Right: Cumulative dose response (cDR) at time T; now the vertical axis shows the integral (area under the curve) z(t)=∫0Ty(t)dt of the response, again plotted against the input u. For example, z1(t) is the area shaded in blue in figure 1. This particular cDR is monotonic. For example, u1<u2<u3 and z1(T)<z2(T)<z3(t), because in figure 1 the area under the red curve y2(t) is larger than the area under the blue curve y1(t), even though the final value y2(T) is smaller than y1(T).

If the turnover of $[eqn]$ is slow, the molecules or other objects represented by $[eqn]$ may accumulate, for instance, in a particular tissue or the bloodstream. It is often the case that one can only measure experimentally, and that a phenotypical response only depends on, the accumulated value ‘or integral under the curve’ (AUC) $[eqn]$ , which we call the cumulative dose response and denote by $[eqn]$ . The right panel of figure 2 shows several cDRs obtained from simulation data as in figure 1. Section 1.9 briefly reviews several areas of biology where cDR’s appear naturally. Specifically, however, this article was motivated by our previous research in immunology that measured cDR’s [6]. We discuss that motivation next.

Motivation for this work: T-cell recognition

1.2.

Adaptation is central to immunity. In particular, T cells must react to stimulation by pathogens and cancers, yet limit their response to maintain self-tolerance and avoid autoimmune reactions. T-cell activation is triggered by the binding to T-cell receptors (TCRs) to peptide major-histocompatibility complex (pMHC) antigens. Activation results in the production of signalling molecules (cytokines), which in turn may recruit other immune components.

The study in [6] examined the response of immune CD8^+^ T cells to external antigen inputs, demonstrating perfect adaptation across a wide range of antigen affinities. Specifically, the experiments in [6] involved stimulating primary human CD8^+^ T cells (with the c58c61 TCR) with recombinant pMHC antigen 9V immobilized on plates, which served as the input $[eqn]$ at various constant concentrations. This antigen is a cancer peptide routinely used in studies of T-cell binding and antigen discrimination. The cumulative amount of cytokine TNF- $[eqn]$ (the output $[eqn]$ ) secreted into the culture medium was measured. Figure 3 shows the cDR when one averages the results of three biological replicates. It shows the average cumulative TNF- $[eqn]$ abundance ( $[eqn]$ ) plotted against several constant concentrations $[eqn]$ , measured at various times ( $[eqn]$ to $[eqn]$ h). Note the non-monotonic, and even somewhat oscillatory, behaviour.

Cumulative dose response based on average of three experiments. Plot uses experimental data from [6] (see also fig. 2B in that paper). Horizontal axis denotes concentrations of the input (in units of ligand in ng/well).

The work in [6] thus raises the question of what network motifs are capable, at least for suitable parameters, to exhibit perfect adaptation as well as non-monotonic cDR as seen in these experiments. It was speculated, on the basis of numerical exploration, that incoherent feedforward loops (IFFLs) cannot result in non-monotonic cDR and thus cannot explain T-cell adaptation as measured by accumulated cytokines, unless a thresholding mechanism is imposed.

Our main results in this article confirm in a mathematically rigorous way that, indeed, the main two common types feedforward loops (called ‘IFFL1’ and ‘IFFL2’ below) can never exhibit such behaviours, because their cDR’s are always monotonic. This is especially surprising for one of them (‘IFFL1’) because for such systems the DR itself can be non-monotonic, yet the cDR is monotonic, as in the cartoon illustrations in figure 2. To see this with an example, consider the following system of two differential equations, which is a particular example of the general equations (1.5) and (1.6) discussed later,

$$[eqn]$$

This is an adapting system: for any given constant input $[eqn]$ , the steady states are $[eqn]$ and $[eqn]$ , which is independent of $[eqn]$ , so that the steady-state output $[eqn]$ is independent of $[eqn]$ . Figure 4 shows plots of the DR (non-monotonic) and the cDR (monotonic).

Plots of DR (y(t)) and cDR (∫0Ty(t)dt) for the example in equations (1.1) and (1.2) . The initial conditions are x(0)=0, y(0)=1/10, and the time horizon is T=1.5. Using logarithmic scale on inputs, for comparison with experimental plots. Observe that the DR is non-monotonic, yet, surprisingly, the cDR is monotonic. Our main theorem proves monotonicity of the cDR in general, for all IFFL1 and IFFL2 systems.

We complement the result for feedforward loops with the new finding that, on the other hand, the standard nonlinear integral feedback for adaptation (‘IFB’ below) is indeed capable of showing non-monotonic cDR, and thus is potentially a mechanism that is consistent with the experimentally observed non-monotonic cDR. To see this with an example, consider the following system of two differential equations, which is a particular example of the general equations (1.9) and (1.10) discussed later,

$$[eqn]$$

This is also an adapting system: for any given constant input $[eqn]$ , the steady states are $[eqn]$ and $[eqn]$ , so that the steady-state output $[eqn]$ is independent of $[eqn]$ . Figure 5 shows plots of the (non-monotonic) cDR.

Plot of cDR (∫0Ty(t)dt) for the IFB example in equations (1.3) and (1.4). The initial conditions are x(0)=0.1, y(0)=6, and time horizons shown are T=3,4,5,6. Observe that, just as with the experimental data plotted in figure 3, the cDR is more monotonic (on the shown ranges, at least) for smaller time horizons T. Using logarithmic scale on inputs, for comparison with experimental plots.

Back to the experimental data, one may ask whether the T-cell experiments point to adaptation, that is, if $[eqn]$ is independent of $[eqn]$ for large $[eqn]$ . Since only $[eqn]$ is experimentally available, we need to estimate the output $[eqn]$ by taking time-derivatives of $[eqn]$ . To obtain a more meaningful estimate than would be obtained from the averages shown in figure 3, we consider instead the separate plots of $[eqn]$ from each experiment. The top panel in figure 6 shows again the cDRs for various times ( $[eqn]$ to $[eqn]$ hours), but now with separate plots from each experiment, starting from the data that were used to generate fig. 1 in the SI of [6].

Top: cDR plots of individual experiments, and measured at different times. Bottom: adaptation behaviour in individual experiments. Output y(t) is estimated from individual cumulative z(t) plots in respective top panels.

Using first-order differences, and imposing a zero value at the start of the experiment: $[eqn]$ , we then derived estimates of $[eqn]$ for the various input values and the three experimental replicates shown in figure 6. See the bottom panel of figure 6. These estimates are very rough, because the time steps are large (1 hour), and in any event experimental data is subject to noise. Nonetheless, these data strongly suggests that adaptation, in the sense that the output $[eqn]$ approaches as time increases a value (in this case $[eqn]$ ), which is the same no matter what the input (drug dose). The estimated negative values of $[eqn]$ at certain time points are likely due to numerical errors or to the fact that there is some cytokine present in the experimental wells which does not arise from the stimulation, so that the values shown are relative to this baseline. In addition to adaptation, the data strongly suggests that the system, or at least its integrated output, is scale-invariant (performs ‘fold change detection’) in the sense of [7,8], at least for large enough input values: the transient outputs are roughly the same (for inputs 3−8), which is a property of systems IFFL2 and the IFB system as discussed below.

Review: network motifs for adaptation

1.3.

As discussed (‘perfect’) adaptation to constant inputs means that, no matter the value of an input, the value of $[eqn]$ for large $[eqn]$ will be the same. Generally speaking, adaptation requires one of two mechanisms for adaptation: incoherent feedforward or negative feedback [9–11].

IFFLs are a type of network motifs that are capable of adaptation [12,13]. In an IFFL, the input $[eqn]$ induces formation of the reporter $[eqn]$ but also acts as a delayed inhibitor, through one or more intermediary control variables. Feedforward motifs are statistically overrepresented in biological systems from bacterial to mammalian cells [14,15]. IFFL’s have been argued to underlie mechanisms involved in such varied contexts as microRNA-mediated loops [16], MAPK pathways [17,18], insulin release [19,20], intracellular calcium release [21,22], Dictyostelium and neutrophil chemotaxis [23,24], NF- $[eqn]$ B activation [25], and microRNA regulation [26], as well as metabolic regulation of bacterial carbohydrate uptake and other substrates [27,28]. IFFL’s may also play a role in immunology, enabling the recognition of dynamic changes in antigen presentation [29] and have been used in synthetic biology to control protein expression under DNA copy variability [30,31]. The paper [32] provided a large number of additional references, and carried out a computational exploration of IFFLs that lead to non-monotonic DRs and/or adaptation.

In IFB loops, the intermediate variable or variables provide a type of memory that integrates the ‘error’ between $[eqn]$ and a steady-state value $[eqn]$ . IFB’s arise in biological systems ranging from E. coli chemotaxis [33] and regulation of tryptophan [34] to human physiological control such as blood calcium homeostasis [35] or neuronal control of the prefrontal cortex [36] to synthetic circuits for adaptation [37–39]. We remark that IFB is in a sense universal for adaptation, because nonlinear changes of coordinates can recast IFFLs that adapt into integral feedback form [40], but these coordinate changes may have no physical interpretation and hence lack interpretability. Moreover, IFBs are known to provide extra degrees of robustness to the adaptation property because, unlike the IFFL, the underlying mechanism can sense and correct for perturbations to the output variable $[eqn]$ . This theme is explored in greater detail in recent papers, where it is shown that IFBs arise naturally when one considers adapting circuits that exhibit a maximal form of robustness [41] or robustness that is independent of the reaction kinetics [42]. The results in these papers apply to arbitrarily-sized networks, and they explicitly identify the network species that create the IFB required for adaptation.

Three paradigmatic circuits with two variables, two types of incoherent loops and one type of IFB, have been much studied, in particular in the context of the ‘fold change’ or ‘scale-invariance’ property in biology [7,8,43], and these systems, discussed next, are the focus of our study (Observe that nonlinearity is essential for non-monotonic cDR’s to be possible, since for linear systems, such as the IFFL $[eqn]$ , $[eqn]$ or the IFB $[eqn]$ , $[eqn]$ , DRs and cDRs will be linear in $[eqn]$ .) Although the concepts introduced here are broadly applicable, we focus our work in these three motifs, but we expect that our results will encourage further work in the direction of characterizing cDR properties for more complex systems.

Three paradigmatic examples of adapting systems

1.4.

In this article, we will consider three types of two-species systems, shown schematically in figure 7. These are the systems studied in [7,8]. In all these examples, $[eqn]$ refers to a positive constant input, $[eqn]$ is the concentration of a ‘controller’ species and the output variable $[eqn]$ is the concentration of the ‘regulated’ or output species.

Three examples of systems: (a) is an IFFL with degradation enhancement, (b) is an IFFL with production inhibition and (d) is an IFB system. The input u is assumed to be a positive constant, and x, y are abundances of a quantity of interest, such as concentration of a protein or mRNA. Arrows ‘→’ indicate positive effect (activation) and blunt edges ‘⊣’ denote negative effects (inhibition). Both (a) and (b) have the same qualitative relation, schematically represented by diagram in (c), between activation of x and y by the input u, and inhibition of y by x, but their dynamical properties are very different. Dynamics can be described by pairs of differential equations for the abundances x(t) and y(t) as a function of time, as shown below the diagrams. State variables x(t) and y(t) are taken to be non-negative. When K=0, it is assumed that x(t)>0 for all t. In these equations, α, β, γ, δ and K are all positive constants.

Steady states and perfect adaptation

1.5.

Let us compute the steady states, obtained by setting the right-hand sides of the differential equations to zero, for the systems shown in figure 7. For the IFFL1 system (a) we have:

$$[eqn]$$ $$[eqn]$$

Finally, for the IFB system (c) we have that there are two types of steady states:

$$[eqn]$$

and, for non-zero $[eqn]$ and assuming $[eqn]$ (for example, if $[eqn]$ ):

$$[eqn]$$

In all three cases, if $[eqn]$ , the system is perfectly adapting.

Scale invariance

1.6.

In addition, when $[eqn]$ both systems IFFL2 and IFB have the scale-invariance (or ‘fold-change detection’, FCD) property [7,8]. This means that the output variable $[eqn]$ satisfies the same differential equation, independently of rescalings of the input $[eqn]$ by any constant factor $[eqn]$ , as shown by the following simple calculation:

$$[eqn]$$

System IFFL1, in contrast, admits no such symmetries.

Stability

1.7.

For both systems IFFL1 (a) and IFFL2 (b) in figure 7, the respective steady states $[eqn]$ are globally asymptotically stable with respect to initial conditions in the positive quadrant $[eqn]$ , $[eqn]$ . This is very simple to show. The variable $[eqn]$ is the solution of a one-dimensional stable linear system, hence converges exponentially to $[eqn]$ . The variable $[eqn]$ is the solution of a time-dependent linear system, with a constant negative rate $[eqn]$ for IFFL2, and a rate for IFFL1 which converges to the strictly negative number $[eqn]$ , and hence also exponentially converges to its steady-state value (see [13] for details, as well as similar results for other IFFL configurations).

The proof of stability for the feedback system IFB (d) in figure 7 requires more work. We proceed by extending the proof from [8], which covered only the case $[eqn]$ . We will assume that $[eqn]$ , which holds in particular if $[eqn]$ . We want to global asymptotic stability with respect to initial conditions with $[eqn]$ positive and $[eqn]$ non-negative (or even arbitrary), when $[eqn]$ is a positive constant, for the two-dimensional system IFB evolving on $[eqn]$ with equations $[eqn]$ , $[eqn]$ .

It is convenient to change coordinates $[eqn]$ and $[eqn]$ , so that we reduce to the study of the system in $[eqn]$ with equations (dropping tildes):

$$[eqn]$$

and we wish to prove the global asymptotic stability of the unique steady state

$$[eqn]$$

Introducing $[eqn]$ , $[eqn]$ , and $[eqn]$ , we can write our system as

$$[eqn]$$

with $[eqn]$ and $[eqn]$ . In other words, we have a mass–spring system $[eqn]$ with nonlinear damping $[eqn]$ and nonlinear spring constant $[eqn]$ . This suggests the use of energy as a Lyapunov function. The map $[eqn]$ is strictly increasing, positive for $[eqn]$ and negative otherwise, and similarly, for $[eqn]$ with respect to $[eqn]$ . Let us define

$$[eqn]$$

By definition, $[eqn]$ and $[eqn]$ for all $[eqn]$ . We also have that $[eqn]$ , and $[eqn]$ , (and mixed partial derivatives are zero), $[eqn]$ is a strictly convex and thus a proper (also called radially unbounded or coercive) function, thus a Lyapunov function candidate. The derivative of $[eqn]$ along trajectories is

$$[eqn]$$

for all $[eqn]$ , and if this derivative only vanishes identically along a trajectory, then $[eqn]$ , which in turn implies, when substituted into $[eqn]$ that $[eqn]$ , i.e. that $[eqn]$ . The LaSalle invariance principle (e.g. [44]) then allows us to conclude global asymptotic stability.

Outline of paper

1.8.

Denoting the input-dependent dynamics as $[eqn]$ for $[eqn]$ , we define the DR and the cDR at time $[eqn]$ as

$$[eqn]$$

In each example, our aim is to determine whether the mapping $[eqn]$ is monotonically increasing. If this monotonicity does not hold universally, we seek to identify sufficient conditions under which it does. It is straightforward to note that if the map $[eqn]$ is monotonically increasing, then the same holds for the map $[eqn]$ . Observe that if the system is linear, both $[eqn]$ and $[eqn]$ will be linear (and therefore monotonic) functions of u. Hence, non-monotonicity is exclusively a property of nonlinear systems.

Let us review the specific examples that we consider. The first is the IFFL1 system shown in figure 7a with equations given by

$$[eqn]$$

with initial state $[eqn]$ and $[eqn]$ . The initial state for $[eqn]$ is arbitrary non-negative, while the initial state for $[eqn]$ is its steady-state value which is independent of $[eqn]$ .

The second example is the IFFL2 system (see figure 7b) given by

$$[eqn]$$

where the initial states $[eqn]$ and $[eqn]$ are arbitrary.

As the final example we have the IFB system (see figure 7c) given by

$$[eqn]$$

with initial states $[eqn]$ (arbitrary) and $[eqn]$ , which is the steady state for $[eqn]$ which is independent of $[eqn]$ .

We remark that for all of these examples, $[eqn]$ (and therefore also $[eqn]$ ) is monotonically increasing for small $[eqn]$ . This is because $[eqn]$ for IFFL1, and $[eqn]$ for IFFL2 and IFB. From $[eqn]$ , it follows that $[eqn]$ and $[eqn]$ , respectively, and both are positive. The situation is far less trivial for larger times $[eqn]$ .

The rest of this article is organized as follows. In the rest of this introduction, we review some other areas of application of cDR’s. In §2, we prove that for IFFL1, even though the map $[eqn]$ may not be monotonic, the map $[eqn]$ is always monotonically increasing, irrespective of the values of $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . In §3, we show that the situation is much simpler for IFFL2, in the sense that both $[eqn]$ and $[eqn]$ are monotonically increasing functions of $[eqn]$ , regardless of the choice of $[eqn]$ , initial states, and the model parameters. Finally, in §4, we show that for the IFB system the map $[eqn]$ is not monotonic in general and find a sufficient condition under which monotonicity holds.

Further motivations for the study of cumulative dose responses

1.9.

In pharmacology and biomedical research, measuring the cumulative amount or ‘area under the curve’ (AUC) of the concentration or abundance of a substance (such as an antibody, cytokine, hormone or metabolite) secreted over a time period, which we termed the cDR, is essential for understanding drug efficacy, toxicity and biological responses. To illustrate the potential wide interest of cDR theory, we briefly discuss here a few examples of such measurements.

Cytokine release assays are laboratory techniques that typically measure the cumulative secretion of cytokines, which are small signalling proteins produced by immune cells in response to stimuli such as drugs, pathogens or immune activators. Cytokines regulate immune and inflammatory responses as well as cell growth, differentiation, survival and tissue repair. Examples are pro-inflammatory (IL-6, IL-1 $[eqn]$ , TNF- $[eqn]$ and IFN- $[eqn]$ ) and anti-inflammatory (IL-10 and IL-4) cytokines, and growth factors (GM-CSF and VEGF). Typical techniques used for measuring cumulative release of cytokines are enzyme-linked immunosorbent assay (ELISA) and multiplex immunoassays such as Luminex $[eqn]$ . For example [45], describes means of measuring pro-inflammatory cytokines in the central nervous system as a way to monitor neuroinflammatory responses to trauma, infection and neurodegenerative diseases. That paper describes an in vivo immunosensing device in which an optical fibre is implanted for a period in the brain of a rodent to capture (by binding to a specific antibody) the cumulative release of a specific cytokine within a region of interest; ELISA is then conducted, in order to determine the cumulative amount of cytokine bound to the fibre. As another example, one of the most common adverse events associated with T-cell bispecific antibody therapies (which are themselves an interesting subject for mathematical modelling [46]) is cytokine release syndrome (CRS), whose symptoms include fever, hypotension, respiratory deficiency and possible multi-organ failure. The paper [47] highlighted the use of Luminex $[eqn]$ and AUC plots of cytokine release to evaluate different therapeutic approaches to the mitigation of CRS.

The metabolic clearance of drugs is often assessed by measuring the cumulative amount of a metabolite over time in biological fluids. Typically, after a drug or prodrug is administered, plasma concentrations are measured as a function of time, and the area under the concentration–time curve is computed [48], thus providing critical insights into drug metabolism, pharmacokinetics and hepatic or renal clearance mechanisms. Hence, helping to understand drug efficacy and safety. The paper [49] provided a systematic evaluation of cumulative drug excretion in clinical pharmacokinetics, emphasizing its role in dosing regimens and safety evaluations. It described various measures of drug accumulation, including the AUC on a graph plotting plasma concentration against time.

As one last example, consider the routine A1c blood test (also known as glycated haemoglobin or HbA1c test), which is key to diabetes diagnosis and management. A1c measures average blood glucose levels over a period of about 2−3 months before testing, and one could therefore think of it as being proportional to the integral of glucose concentration over that period. The integration effect, compared with measuring of short-term fluctuations in blood glucose, is due to the long lifespan of red blood cells, whose haemoglobin binds glucose to form HbA1c.

Unconditional monotonicity of the cumulative dose response for incoherent feedforward loop 1

We shall prove the monotonicity of the map $[eqn]$ in multiple steps. As the first step, we simplify the system (1.5) and (1.6) in §2.1 to obtain a parameter-free form that is more amenable to analysis. In §2.2, we then derive explicit expressions for the cDR and its partial derivative with respect to $[eqn]$ . This partial derivative is given by an integral expression, and in §2.3, we establish a couple of properties of the integrand and the asymptotic value of the integral expression. These results enable us to prove the monotonicity of the cDR in §2.3 by demonstrating that its partial derivative with respect to the input $[eqn]$ is always non-negative.

Simplifying the system

2.1.

Let $[eqn]$ be the solution of (1.5) and (1.6) for $[eqn]$ . We scale time by $[eqn]$ and the state values by a suitable ratio to define

$$[eqn]$$

Then, the dynamics of this rescaled system are given by

$$[eqn]$$

Therefore, if we define

$$[eqn]$$

then $[eqn]$ satisfies the ODEs

$$[eqn]$$

with initial states $[eqn]$ and $[eqn]$ . Let $[eqn]$ be the solution of this system. To prove the result it suffices to show that the map

$$[eqn]$$

is monotonically increasing for any $[eqn]$ .

Henceforth, we shall drop the hats for notational convenience and suppose that the dynamics is given by

$$[eqn]$$

with initial state $[eqn]$ (arbitrary) and $[eqn]$ . Letting $[eqn]$ be the solution of this system, we shall show that the cDR map

$$[eqn]$$

is monotonically increasing for any $[eqn]$ . In order to prove this we will prove that for any $[eqn]$ and $[eqn]$

$$[eqn]$$

where $[eqn]$ denotes the partial derivative with respect to $[eqn]$ .

Derivation of explicit expressions

2.2.

We now develop explicit expressions for $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . It is easy to see that

$$[eqn]$$

which also implies that

$$[eqn]$$

Note that ODE (2.2) can be written as

$$[eqn]$$

Multiplying with the integrating factor $[eqn]$ on both sides we obtain

$$[eqn]$$

Integrating both sides and using $[eqn]$ we get the usual variation of the parameters formula:

$$[eqn]$$

From (2.1), we know that

$$[eqn]$$

where the second relation follows from (2.4). Plugging this in the previous expression for $[eqn]$ we get

$$[eqn]$$

where to derive the last expression we have made a change of variable from $[eqn]$ to $[eqn]$ .

Note that by changing the order to integration in the second term in the r.h.s. below we obtain

$$[eqn]$$

Making the change of variable $[eqn]$ we see that

$$[eqn]$$

where $[eqn]$ is the special exponential integral function defined as

$$[eqn]$$

Plugging this in we obtain

$$[eqn]$$

where

$$[eqn]$$

Observe that at $[eqn]$ we have $[eqn]$ and we claim that as $[eqn]$ , we have $[eqn]$ . To see this note that as $[eqn]$ , both $[eqn]$ and $[eqn]$ approach $[eqn]$ and close to 0 we have

$$[eqn]$$

which is just $[eqn]$ as $[eqn]$ . By the definition of the Ei function (2.7)

$$[eqn]$$

and hence by the chain rule

$$[eqn]$$

and

$$[eqn]$$

Using these expressions we can write the derivative of $[eqn]$ as

$$[eqn]$$

Therefore, applying integration by parts to the second term in the r.h.s of (2.8) we get

$$[eqn]$$

Upon substituting this term in the r.h.s. of (2.8) we obtain

$$[eqn]$$

Note that by rearranging and simplifying, (2.9) can be expressed as:

$$[eqn]$$

Differentiating (2.10) with respect to $[eqn]$ we obtain

$$[eqn]$$

which can be rewritten as

$$[eqn]$$

Two useful results

2.3.

We now establish a couple of useful results that will help us in proving the monotonicity property of the cDR.

Lemma 2.1. Fix any positive $[eqn]$ and $[eqn]$ , and for $[eqn]$ define functions

$$[eqn]$$

and

$$[eqn]$$

Then, the function $[eqn]$ can cross the $[eqn]$ -axis only once in the interval $[eqn]$ .

Proof. Note that $[eqn]$ . Setting $[eqn]$ we get

$$[eqn]$$ $$[eqn]$$

Let $[eqn]$ be the function

$$[eqn]$$

To prove the lemma it suffices to show that the function $[eqn]$ is monotonically increasing, which we shall show by proving that $[eqn]$ for all $[eqn]$ . For convenience let us define a function

$$[eqn]$$

Then, we can write $[eqn]$ as

$$[eqn]$$

which shows that

$$[eqn]$$

Suppose we can show that

$$[eqn]$$

Then, for $[eqn]$ we have $[eqn]$ and so $[eqn]$ . On the other hand, for $[eqn]$ , $[eqn]$ and using $[eqn]$ , we obtain

$$[eqn]$$

Hence, to prove the lemma we just need to prove the inequality (2.13) to which we come to now. Note that

$$[eqn]$$

Since $[eqn]$ for any $[eqn]$ , we would have that $[eqn]$ . Since the function $[eqn]$ is monotonically decreasing for $[eqn]$ we must have

$$[eqn]$$

Observe that we can write $[eqn]$ as

$$[eqn]$$

As the first term is always positive, we have $[eqn]$ , which along with (2.14) shows that

$$[eqn]$$

This completes the proof of the inequality (2.13) and concludes the proof of this lemma.∎

Proposition 2.2. For any $[eqn]$ , the integral

$$[eqn]$$

has a positive value.

Proof. Define a function

$$[eqn]$$

Note that the integral (2.15) can be rewritten as

$$[eqn]$$

where

$$[eqn]$$

Hence, in order to prove that $[eqn]$ is positive, we just need to prove that the function $[eqn]$ is monotonically increasing. This is what we show next.

Using the fact that

$$[eqn]$$

we obtain

$$[eqn]$$

Applying the change of variable $[eqn]$ , we get

$$[eqn]$$

The integral on the right can be expressed in terms of the Gamma function $[eqn]$ as

$$[eqn]$$

Substituting this we obtain

$$[eqn]$$

Since $[eqn]$ we can express $[eqn]$ as

$$[eqn]$$

where

$$[eqn]$$

As each $[eqn]$ is a product of positive monotonically increasing functions, the function $[eqn]$ is also monotonically increasing. This completes the proof of this proposition.∎

Proving monotonicity of the cumulative dose response

2.4.

Finally, we now prove the monotonicity of the cDR by proving (2.3). For this we shall use the integral expression (2.11).

Let us first deal with the case $[eqn]$ . In this scenario, $[eqn]$ and since $[eqn]$ we have

$$[eqn]$$

Therefore, the integrand on the r.h.s. of (2.11) can be bounded below

$$[eqn]$$

Hence, the integrand in (2.11) is always non-negative which establishes the monotonicity of cDR for $[eqn]$ .

We now consider the case $[eqn]$ . Note that in this case $[eqn]$ (see (2.4)) is monotonically increasing from $[eqn]$ to $[eqn]$ . Hence, $[eqn]$ and

$$[eqn]$$

which allows us to conclude that $[eqn]$ for all $[eqn]$ , by a simple comparison argument. Suppose that $[eqn]$ . In this case, we prove that for any positive $[eqn]$ , the trajectory $[eqn]$ can only change its sign once, to go from positive to negative, and then it will stay negative and asymptotically approach $[eqn]$ . To see this define

$$[eqn]$$

and then

$$[eqn]$$

Since we have assumed $[eqn]$ we have the differential inequality

$$[eqn]$$

Hence, by the comparison argument, if there exists a $[eqn]$ such that $[eqn]$ , then for all $[eqn]$ we have $[eqn]$ which also implies that

$$[eqn]$$

Now let $[eqn]$ be the first zero of the trajectory $[eqn]$ . Observe that

$$[eqn]$$

As $[eqn]$ is the first zero of the trajectory $[eqn]$ we have $[eqn]$ for $[eqn]$ , and $[eqn]$ and $[eqn]$ . Hence, we must have $[eqn]$ and by previous arguments $[eqn]$ for all $[eqn]$ . Therefore, in this interval $[eqn]$ , we have the inequality

$$[eqn]$$

and since $[eqn]$ , the comparison argument shows that $[eqn]$ for all $[eqn]$ . Hence, the trajectory $[eqn]$ can only change its sign once, to go from positive to negative. Therefore, in order to prove $[eqn]$ for any $[eqn]$ , it suffices to prove that this holds in the limit $[eqn]$ . Letting $[eqn]$ in (2.11) we arrive at

$$[eqn]$$

Since $[eqn]$ for all $[eqn]$ , in order to prove the positivity of this integral it suffices to prove its positivity for $[eqn]$ , i.e.

$$[eqn]$$

which holds due to Proposition 2.2.

We now come to the case where $[eqn]$ . Set $[eqn]$ and let $[eqn]$ and $[eqn]$ be the functions defined in the statement of Lemma 2.1. Note that (2.11) can be written in terms of this function as

$$[eqn]$$

Moreover

$$[eqn]$$

Lemma 2.1 proves that the function $[eqn]$ can only cross the $[eqn]$ -axis at most once in the interval $[eqn]$ . If it does cross then it goes from positive to negative and it stays negative untill it becomes 0 at $[eqn]$ . The same holds for the function $[eqn]$ . Hence, if we can prove that

$$[eqn]$$

then, as $[eqn]$ , it automatically implies that

$$[eqn]$$

To see this let $[eqn]$ be the time in $[eqn]$ where the function $[eqn]$ crosses the $[eqn]$ -axis. Then (2.18) implies that

$$[eqn]$$

Therefore, using that $[eqn]$ for $[eqn]$ and $[eqn]$ for $[eqn]$ , we get

$$[eqn]$$

This shows that to prove the monotonicity result it suffices to prove (2.18), which is of course equivalent to proving the positivity of (2.11) for $[eqn]$ and $[eqn]$ . As mentioned above, this positivity follows from (2.16) which is shown in Proposition 2.2. This completes the proof of the cDR monotonicity result for IFFL1.

Unconditional monotonicity of both dose response and cumulative dose response for incoherent feedforward loop 2

We now consider the IFFL2 system described by equations (1.7) and (1.8) . For this system, we can prove that the DR map, $[eqn]$ , is a monotonically increasing function of the input $[eqn]$ , and hence, the same holds for the cDR map. A direct proof of this is provided in §3.1, while a more conceptual argument, based on the theory of monotone systems, is presented in §3.2. Although the former approach is simpler for this particular example, the latter is more generalizable to other examples.

Direct proof

3.1.

In order to prove the monotonicity of the map $[eqn]$ , it suffices to show that the partial derivative $[eqn]$ is non-negative for any $[eqn]$ and $[eqn]$ . Since $[eqn]$ satisfies the linear ODE (1.7), we can solve for it explicitly to obtain

$$[eqn]$$

which also implies that

$$[eqn]$$

As $[eqn]$ satisfies the ODE (1.8), we can differentiate it with respect to $[eqn]$ to obtain an ODE for $[eqn]$ as

$$[eqn]$$

Substituting $[eqn]$ from (3.1) and re-arranging we get

$$[eqn]$$

Since $[eqn]$ , it is immediate that $[eqn]$ , and so we have the differential inequality

$$[eqn]$$

As $[eqn]$ , by the comparison argument it follows that $[eqn]$ for all $[eqn]$ and $[eqn]$ . This concludes the proof of the monotonicity result for IFFL2.

Monotone systems: an approach to show dose-response monotonicity of incoherent feedforward loop 2

3.2.

Monotone systems were introduced in the pioneering work of Smale, Smith, Hirsch, Mallet-Paret, Sell and others [50–53]. They have the property that larger initial conditions give rise to larger state trajectories, where ‘larger’ is interpreted according to a specified order in the state variables. Special cases of monotonic systems are obtained when the order is a coordinatewise order. For example, in two dimensions, the ‘NorthEast’ (NE) is defined by saying that a point $[eqn]$ is larger than a point $[eqn]$ if both $[eqn]$ and $[eqn]$ , that is, if it is to the ‘north’ and ‘east’ (higher, to the right) of the second point; similarly the ‘NorthWest’ order would be defined by asking that $[eqn]$ and $[eqn]$ . (Note that these are ‘partial orders’ in the sense that two vectors may not be comparable: for example, neither $[eqn]$ nor $[eqn]$ is larger than the other in the NE order.) The generalization to external inputs and outputs [54] enabled the development of a network interconnection theory as well as leading to conclusions regarding the effect of inputs: for example, monotonic inputs result in monotonic transient behaviour [55]. This means that the DR (and thus also the cDR) will always be monotonic, for monotone systems.

These developments led to monotone systems playing a key role in analysing the global behaviour of dynamical systems in various areas of engineering as well as biology [56]. What makes monotone system theory so useful is that there are ways to check monotonicity without solving a set of differential equations $[eqn]$ . For example, for the $[eqn]$ -dimensional analogue of the NE order one requires that the off-diagonal terms of the Jacobian matrix of $[eqn]$ should all be non-negative, and a similar condition holds if there are inputs. More generally, monotonicity with respect to some (not necessarily the NE) coordinatewise order requires that all loops in the interaction graph obtained from the off-diagonal terms of the Jacobian should have a net positive sign (see [57] for more discussion and examples).

When trying to apply monotone systems theory to our three paradigmatic examples IFFL1, IFFL2 and IFB, an obvious problem arises: these systems are not monotone with respect to any possible coordinatewise order, as IFFLs and negative feedback loops contradict the positive loop condition. However, in the case of the IFFL2 system, which we reproduce here for convenience:

$$[eqn]$$

there is a Lie group of symmetries or ‘equivariances’ that preserve the output. These equivariances were the main object of study in the work in [8] on scale invariance, and were key to the analysis of an immunology model in [29]. Specifically, when $[eqn]$ the discussion in §1.6 showed that the equations do not change under the one-parameter Lie group of transformations $[eqn]$ , and in particular scaling $[eqn]$ and $[eqn]$ by the same constant does not alter the dynamics of $[eqn]$ . This suggests introducing the new variable $[eqn]$ . Using the variables $[eqn]$ and $[eqn]$ , the equations become:

$$[eqn]$$

This is a monotone system, because the only off-diagonal term in the Jacobian is $[eqn]$ . Therefore, the trajectories depend monotonically on $[eqn]$ , and hence also depend monotonically on $[eqn]$ . A similar argument works for $[eqn]$ , but now the FCD property fails and the equivariance will provide merely an embedding into a monotone system rather than an equivalence. Indeed, let $[eqn]$ . Now the $[eqn]$ equation is no longer decoupled from $[eqn]$ . However, we can look at the following extended system (we add an equation $[eqn]$ to convert the external input into a state variable):

$$[eqn]$$

The off-diagonal elements of the Jacobian are $[eqn]$ , $[eqn]$ , and $[eqn]$ , all positive. Thus, the extended system is monotone, and therefore all variables, and in particular $[eqn]$ , depend monotonically on $[eqn]$ . This shows monotonicity of the DR, as claimed.

Conditional monotonicity of the cumulative dose response integral feedback

As our final example, we examine the IFB system defined by equations (1.9) and (1.10). Our objective is to show that, unlike the previously discussed incoherent feedforward loop systems, this negative feedback system can exhibit non-monotonic behaviour in the mapping $[eqn]$ . We also identify a condition under which this mapping remains monotonic. To simplify the analysis, we restrict our attention to the case $[eqn]$ . If the property of monotonicity or non-monotonicity arises in this case, it will also persist for small positive values of $[eqn]$ due to the parametric continuity of solutions.

We begin our analysis by simplifying the equations for the IFB system using scaling arguments in §4.1. Here, we also reformulate the dynamics in terms of the ‘error’, which is obtained by subtracting the input-independent steady state from $[eqn]$ . Notably, the monotonicity properties of the original cDR map can be equivalently studied using the error cDR map. In §4.2, we prove that the partial derivative of the error cDR map, with respect to the input value $[eqn]$ , is always positive at steady state, confirming that the steady-state cDR map is monotonic. For finite-time analysis, we establish a connection between the dynamics of the error and a ‘damped’ harmonic oscillator in §4.3. This connection allows us to derive a conditional monotonicity result for the error cDR map in §4.4. Furthermore, this result enables us to identify a finite range of time values for each input value $[eqn]$ , within which the error cDR map may become non-monotonic. We leverage this insight to numerically illustrate this non-monotonicity in §4.5.

Simplifying the system

4.1.

Recall the IFB equations (1.9) and (1.10) . We scale time by $[eqn]$ and the state values by a ratio to define

$$[eqn]$$

Then, we can write this system as

$$[eqn]$$

Therefore, if we define

$$[eqn]$$

then $[eqn]$ satisfies the ODEs

$$[eqn]$$

with initial values $[eqn]$ and $[eqn]$ . Henceforth, we shall drop the hats for notational convenience and suppose that the dynamics is given by

$$[eqn]$$

with initial values $[eqn]$ and $[eqn]$ . Let $[eqn]$ be the solution of this initial value problem for $[eqn]$ .

Note that $[eqn]$ is both the initial state and the steady-state value of the output variable $[eqn]$ . We shall reformulate the dynamics in terms of the ‘error’

$$[eqn]$$

and variable $[eqn]$ defined as

$$[eqn]$$

Then, dynamics of $[eqn]$ is given by

$$[eqn]$$

Hence, we can solve for $[eqn]$ as

$$[eqn]$$

and the dynamics of the error is given by

$$[eqn]$$

Henceforth, instead of examining the monotonicity of the original cDR map, we shall equivalently examine the monotonicity of the error cDR map given by

$$[eqn]$$

Steady-state analysis

4.2.

We proved earlier that the steady state is stable. Thus, $[eqn]$ and so $[eqn]$ . Therefore, expression (4.3) implies that

$$[eqn]$$

which is a monotonically increasing function of $[eqn]$ . In particular, by differentiating with respect to $[eqn]$ we obtain

$$[eqn]$$

If this positivity holds for any finite time-interval then we shall have the monotonicity of the error cDR map. We shall show that this monotonicity does not always hold and identify a sufficient condition under which it does. For this we shall connect this problem to a harmonic oscillator with a time-varying frequency.

Connection to a harmonic oscillator

4.3.

Differentiating the error equation (4.4) with respect to $[eqn]$ and using that $[eqn]$ we get

$$[eqn]$$

which means that the error $[eqn]$ satisfies the equation for a ‘damped’ harmonic oscillator with a non-constant frequency:

$$[eqn]$$

with initial conditions $[eqn]$ and $[eqn]$ . Let us define the ‘frequency’ and Hamiltonian for this oscillator as

$$[eqn]$$

respectively. Then, their dynamics can be derived as

$$[eqn]$$

Using (4.7) we obtain

$$[eqn]$$

Observe that the error $[eqn]$ in our adapting circuit goes to $[eqn]$ as $[eqn]$ . Equation (4.9) shows that when this error is below $[eqn]$ , the Hamiltonian is decreasing. Note that the initial value of this Hamiltonian is

$$[eqn]$$

This brings us to a proposition that shows that if this value is less than $[eqn]$ , then the error $[eqn]$ remains below $[eqn]$ at all times.

Proposition 4.1 Suppose that the following holds:

$$[eqn]$$

Then we must have that $[eqn]$ for all $[eqn]$ .

Proof. Condition (4.10) implies that $[eqn]$ . Recall that $[eqn]$ . Let $[eqn]$ be the first time $[eqn]$ reaches $[eqn]$ , i.e. $[eqn]$ and $[eqn]$ for all $[eqn]$ . Then due to (4.9) the Hamiltonian must be decreasing in the interval $[eqn]$ , and so $[eqn]$ . But this is a contradiction since $[eqn]$ and so by definition $[eqn]$ . Hence, $[eqn]$ which means that $[eqn]$ for all $[eqn]$ .∎

Conditional monotonicity result

4.4.

Now let us substitute $[eqn]$ from (4.3) into equation (4.4) to obtain

$$[eqn]$$

Differentiating this with respect to $[eqn]$ we get

$$[eqn]$$

This shows that if we define $[eqn]$ as

$$[eqn]$$

then $[eqn]$ also satisfies the following second-order ODE for the damped harmonic oscillator (4.7) with initial conditions $[eqn]$ and $[eqn]$ . To prove the monotonicity we need to show that

$$[eqn]$$

which is equivalent to proving that

$$[eqn]$$

Again we can definite the Hamiltonian as

$$[eqn]$$

and it will have the same dynamics as before:

$$[eqn]$$

which we can also write as

$$[eqn]$$

We now come to our main result for the IFB example, which proves monotonicity of the error cDR under condition (4.10).

Proposition 4.2 Suppose that condition (4.10) holds. Then, the map

$$[eqn]$$

is monotonically increasing for any $[eqn]$ .

Proof. When (4.10) holds, we have that the error $[eqn]$ for all $[eqn]$ due to proposition 4.1. This fact along with (4.14) implies that $[eqn]$ for $[eqn]$ . Therefore, (4.12) holds which proves the required monotonicity property.∎

Numerical illustration of non-monotonicity

4.5.

When condition (4.10) does not hold, we do not have any analytical approach for checking the monotonicity property. Therefore, we need to rely on numerical simulations. For any given $[eqn]$ , the monotonicity condition (4.12) can fail at any $[eqn]$ . However, as we cannot check this condition for all $[eqn]$ , we first show that if this monotonicity fails then it has to fail in the finite time - interval $[eqn]$ where $[eqn]$ is defined by

$$[eqn]$$

Note that $[eqn]$ is finite because both $[eqn]$ and $[eqn]$ are non-negative quantities converging to zero as $[eqn]$ . This is due to stability and the fact that (4.6) implies that $[eqn]$ as $[eqn]$ .

Proposition 4. $[eqn]$ be the finite time defined by (4.16). Then, we must have $[eqn]$ for all $[eqn]$ which is the same as saying that $[eqn]$ for all $[eqn]$ .

Proof. Pick a small $[eqn]$ and define

$$[eqn]$$

Note that for the same reason that $[eqn]$ is finite, $[eqn]$ must also be finite. Moreover, as $[eqn]$ decreases the set $[eqn]$ gets bigger and consequently its infimum, which is $[eqn]$ , gets smaller. Therefore, $[eqn]$ decreases monotonically as $[eqn]$ decreases and since $[eqn]$ and $[eqn]$ are continuous functions of time we must have

$$[eqn]$$

By the definition of $[eqn]$ we have that $[eqn]$ which implies that $[eqn]$ . Let $[eqn]$ be the first time after $[eqn]$ such that $[eqn]$ and consequently $[eqn]$ . On the open interval $[eqn]$ , $[eqn]$ must be decreasing due to (4.9) and so we should have $[eqn]$ . However, this leads to a contradiction because $[eqn]$ . Therefore, $[eqn]$ which means that $[eqn]$ for all $[eqn]$ . This also means that $[eqn]$ is decreasing in the interval $[eqn]$ due to (4.14). Since $[eqn]$ this implies that $[eqn]$ for all $[eqn]$ . Using the continuity of $[eqn]$ and the limit (4.17) we can conclude that $[eqn]$ for all $[eqn]$ . This completes the proof of this proposition.∎

Figure 8 illustrates this proposition for three values of $[eqn]$ . One can observe that $[eqn]$ holds for all $[eqn]$ for each $[eqn]$ -value.

Plot of the map t↦u∫0t∂uy(s)ds for three values of u. Note that for each u, we have u∫0t∂uy(s)ds≥0 for all t≥Tu.

We can now numerically solve the system $[eqn]$ in the interval $[eqn]$ and check if $[eqn]$ exceeds $[eqn]$ or not. For any $[eqn]$ such that $[eqn]$ we would have $[eqn]$ , which would imply non-monotonicity. For any given $[eqn]$ we define the non-monotonicity score as

$$[eqn]$$

where $[eqn]$ (resp. $[eqn]$ ) denotes the positive (resp. negative) part of $[eqn]$ . In figure 9 we plot the non-monotonicity score $[eqn]$ and the time $[eqn]$ for a range of $[eqn]$ -values for $[eqn]$ . Notice that $[eqn]$ is monotonically increasing with $[eqn]$ , while the non-monotonicity score $[eqn]$ stays at zero untill $[eqn]$ and then it gradually starts increasing. This shows that for higher values of $[eqn]$ the map

Plot of the maps u↦S(u) and u↦T(u) with x0=p=1.0.

$$[eqn]$$

is always non-monotonic for some $[eqn]$ which will increase as $[eqn]$ increases. Figure 10 illustrates this non-monotonicity for some values of $[eqn]$ and $[eqn]$ .

Plot of the map u↦∫0Ty~u(t)dt for some values of p and T. Note that this map exhibits non-monotonicity in these cases.

Conclusions

We introduced the notion of cDR, and went on to show mathematically that both the incoherent feedback loop IFFL1 and IFFL2 motifs can only produce monotonic cDRs, even if (for IFFL1) the DR itself may be non-monotonic. On the other hand, we also established conditions under which the IFB mechanism can produce a non-monotonic cDR.

Leaving aside linear systems, which can never lead to non-monotonic (or, for that matter, any nonlinear) cDR, the motifs that we analysed are considered the simplest paradigms for adaptation in biology [7,8,43]. The concepts introduced here are broadly applicable, even if results were established for idealized two-variable systems that adapt perfectly. They provide a foundation for developing a more comprehensive mathematical theory that can qualitatively characterize cDR maps in more complex, multi-variable systems. One particularly interesting direction would be to study how network interconnections (such as cascades) of these motifs can preserve the qualitative properties of cDR’s. Another is to include scenarios with non-ideal behaviours, such as species dilution, saturation in reaction rates and resource competition, as well as imperfect adaptation.

One may view our study of the cDR properties as an addition to the toolkit of mathematical methods for model discrimination and invalidation, in a spirit similar to the work in [58] that infers the existence of IFFLs or negative feedback loops when time responses are non-monotonic, or to the work in [10] on using periodic inputs to rule out IFFLs as the basis of adaptation. In such a role, one can perform dose-dependent experiments and on the basis of plotted cDR’s, eliminate systems whose structure is not consistent with an observed non-monotonic cDR. Conversely, one could ask how to introduce new variables or otherwise modify a model to match such qualitative behaviours.