Linear models of activation cascades: analytical solutions and coarse-graining of delayed signal transduction

Cellular signal transduction usually involves activation cascades, the sequential activation of a series of proteins following the reception of an input signal. Here, we study the classic model of weakly activated cascades and obtain analytical solutions for a variety of inputs. We show that in the special but important case of optimal gain cascades (i.e. when the deactivation rates are identical) the downstream output of the cascade can be represented exactly as a lumped nonlinear module containing an incomplete gamma function with real parameters that depend on the rates and length of the cascade, as well as parameters of the input signal. The expressions obtained can be applied to the non-identical case when the deactivation rates are random to capture the variability in the cascade outputs. We also show that cascades can be rearranged so that blocks with similar rates can be lumped and represented through our nonlinear modules. Our results can be used both to represent cascades in computational models of differential equations and to fit data efficiently, by reducing the number of equations and parameters involved. In particular, the length of the cascade appears as a real-valued parameter and can thus be fitted in the same manner as Hill coefficients. Finally, we show how the obtained nonlinear modules can be used instead of delay differential equations to model delays in signal transduction.


Introduction
Activation cascades are pervasive in cellular signal transduction systems [1,2]. In its simplest form, an activation cascade comprises a set of components (typically proteins) that become sequentially activated in response to an external stimulus (figure 1). These systems have been the subject of numerous studies, experimental and theoretical [1,[3][4][5][6][7][8][9]. The role of activation cascades in cellular signal transduction is manifold. Cascades can relay, amplify, dampen or modulate signals in order to achieve a variety of cellular responses. One of the best-studied examples of such a system is the mitogen-activated protein kinase (MAPK) cascade, which plays a central role in key cellular functions, such as regulation of the cell cycle, stress responses and apoptosis [2].
Models of activation cascades are known to exhibit a range of nonlinear behaviours, including ultrasensitivity [6,10] and multistability [5,11]. Linearized models of cascades [1] (the so-called 'weakly activated' regime studied here) are also of theoretical interest, and have been studied to evaluate signalling times [12], signal specificity [13] and optimal gain [4]. Such linearized descriptions of cascades often appear as part of larger and more complicated models, and have been shown to be useful in model reduction techniques [14]. Hence obtaining coarse-grained representations of such cascades would be useful not only to simplify their mathematical analysis but also computationally, to allow for compact implementations in models for systems biology. Furthermore, weakly activated cascades are of importance in quantitative biology as they have been observed experimentally [15]. In this context, it would be desirable to estimate the length of an unobserved cascade from data without having to create and fit several models, each with a different number of equations to represent the varying length of the cascade.
Here, we present a study of analytical solutions of ordinary differential equation (ODE) models of linear activation cascades. First, we obtain general solutions for weakly activated cascades. We then focus on the case when the gain of the cascade is optimal (i.e. when all deactivation rates are identical), and find that a lower incomplete gamma function with only three real-valued parameters represents the output of the entire cascade. We exemplify the use of this coarsegrained solution to describe the downstream output induced by several time-dependent inputs of interest, including step functions, exponentially decaying signals, Gaussian inputs and periodic stimuli. We also show that the obtained solution has real-valued parameters directly linked to the length and filtering properties of the cascade, and can thus be used to fit data capturing efficiently the delay and distortion introduced by the cascade. We also explore the application of our results to non-optimal cascades, i.e. when the requirement of identical deactivation rates is relaxed. When only one deactivation rate is different, the equations can be reordered, so that a lumped gamma function representation can be used for the block of identical proteins without altering the final output of the cascade. We also show that when the deactivation rates are randomly distributed, the gamma function can still be used to represent the distribution of the outputs of the cascade. Finally, we show how the gamma function representation of a cascade can be used as a computationally efficient replacement of delay differential equations (DDEs).

Weakly activated cascades and their gamma function solution
Consider a cascade involving n components that are activated in succession. Upon perception of the input signalR(t), the first inactive component (x 1 ) is transformed into its activated form (x Ã 1 ), which then activates the next component (x 2 ). Sequential activation of x i by x Ã iÀ1 continues until the end of the cascade. The output of the cascade of length n is the activated form of the last component, x Ã n . In the case of the MAPK cascade, the components are proteins, and the activation corresponds to a post-translational modification, i.e. phosphorylation. However, the formalism can also describe other sequential biochemical processes with similar functional relationships, e.g. n-step deoxyribonucleic acid (DNA) unwinding [16].
If we use mass-action kinetics without an intermediate complex to describe protein activation, the reaction describing the activation of x 1 is and for the rest of the proteins x i (i ¼ 2, . . . ,n) we have We also assume that all proteins deactivate spontaneously with constant rate The system of nonlinear ODEs describing the time evolution of the full activation cascade is [1] dx Ã where we have defined the total amount of each protein We also assume that the model operates over time scales where there is no significant protein production, so that the amount of each protein T i can be considered constant. If the time scales are such that the total amount of protein varies significantly, then each T i would have to be described by its own ODE according to additional biological knowledge.

The general solution for weakly activated cascades
As shown in [1], in the weakly activated regime T i ) x Ã i , one takes the approximation (T i À x Ã i ) % T i , and the original system (2.1) can be rewritten as a driven linear system: 2Þ . , 0 T is the first nÂ1 vector of the canonical basis, and the n Â n rate matrix A is . . . a n Àb n where a i ¼â i T i , 8i. This system can be solved using the Laplace transform with auxiliary variable s. If the cascade receives an integrable input R(t), it is easy to show that the Laplace transform of the kth protein is

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{
Correction for initial condition : ð2:4Þ The first term on the right-hand side corresponds to the Laplace transform of x Ã k (t) for initial conditions the cascade starts from rest), and the input: Figure 1. A typical protein activation cascade of length n. The proteins (nodes) in the cascade can either be in an inactive (x i ) or active (x Ã i ) state. An external signal R(t) activates the first node. Once a node is active, it activates the next component in the cascade until the end. The activation rate of each x i is a i , and the deactivation rate of each x Ã i is b i . Adapted from [1]. (Online version in colour.) rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160409 second term contains the correction for non-zero initial conditions. The term a (k) is the geometric mean of the activation rates up to k: : ð2:6Þ Using linearity and the convolution properties of the Laplace transform, the output of the cascade is finally obtained as and the pre-factor incorporates the product of all the activation rates, Although equation (2.7) describes the evolution of a general initial condition, in this study we will assume henceforth that the cascade is initially fully inactive (i.e. x Ã i (0) ¼ 0, 8i). In the cases when x Ã i (0) = 0, then the exponential correction introduced by the initial conditions can be incorporated to the calculations.
Example 2.1. If a linear cascade is subject to a constant stimulus given by the step function R(t) ¼ 1, t ! 0, and shows that the output of the last protein in the cascade is given by ð2:8Þ

Optimal linear cascades
Activation cascades are substantial modules of the cellsignalling machinery and, as such, they should be efficient in minimizing the use of energetic resources, such as adenosine triphosphate, or of cellular building blocks, such as amino acids. In the study of Chaves et al. [4], it was shown that when a weakly activated cascade (2.2) is required to provide a given gain, the amplification is achieved optimally when the number of steps in the cascade (e.g. the number of proteins) is finite and all deactivation rates are equal, i.e. b i ¼ b, 8i. This result means that arbitrarily long cascades are not useful for cells when a particular amplification gain from external signals is required. For an optimal cascade, the rate matrix in equation (

Linear cascades under different input functions
We now consider the time-dependent output of a cascade under four different inputs of biological interest.

Step-function stimulus
In an experimental setting, one often studies the response of a biological system to a step-function stimulus such as constant temperature, light or treatment started at time t ¼ 0. In this case, the stimulus is where I n is the n Â n identity matrix, and e tA is the matrix exponential.
If the cascade is optimal (i.e. A ¼Ã), the Laplace transform of the last protein given by (2.4) becomes and taking the inverse transform we obtain is the normalized lower incomplete gamma function whose general form is [17] where G(a) is the gamma function and e Às s aÀ1 ds, Re(a) . 0:

Exponentially decreasing stimulus
When the first protein in the cascade is subject to an exponentially decaying stimulus (e.g. when the input is a reactive molecule or a molecule that becomes metabolized, or if the receptors become desensitized) rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160409 If we assume that the cascade is optimal (A ¼Ã), then and the output of the cascade is given by where a (n) is defined in (2.5). As in the case of constant stimulus, the solution is also given in terms of the lower incomplete gamma function.

Periodic stimulus
In certain experimental settings, we are interested in the response of a system to a periodic stimulus, e.g. circadian rhythms or day/night cycles [18]. Let us consider a linear cascade of length n with periodic input which oscillates between 0 and 2 with mean 1 and frequency v . 0. From a resting initial condition, the solution to equation (2.2) is When the cascade is optimal (A ¼Ã), the explicit solution for the nth protein in the cascade is , u ¼ arctan (b=v) and the T nþk ( cos u) are the Chebyshev polynomials evaluated at cosu.
Asymptotic limits provide useful insights. When the frequency is large compared with the deactivation rate, i.e. v ) b, then b=r ≃ 0, u ≃ 0 and we obtain Hence for large frequencies the oscillations in equation (3.8) are filtered out, and the solution approaches the response to the step function given by equation (3.2). Conversely, when the deactivation of the proteins dominates the frequency (i.e. b ) v) the behaviour of x Ã n will be dominated by the sinusoidal input.
In general, asymptotically as t ! 1, the cascade acts broadly as a filter with an overall amplification (a (n) =b) n , and an oscillatory term attenuated by a factor (b=r) n with a delay phase nu: where we have used the fact that lim t!1 P(n, bt) ¼ 1. Note . 0 for all t. Cascades with more complicated temporal stimuli can be analysed similarly using the Fourier series expansion of R(t).

Gaussian stimulus
Gaussian input functions are employed to represent drug intake and other such signals. Consider a cascade of length n with input which describes a bell-curve centred at t ¼ m, with height 1 and amplitude z. The solution of equation (2.2) from inactive initial conditions under this input is then given by When a Gaussian input (2ps 2 ) À1=2 e À((tÀm) 2 =2s 2 ) becomes increasingly narrow (i.e. s ! 0), it approaches in the limit a Dirac delta function: In that case, from equation (2.7) the solution for the nth protein is 4. Applications of the analytical solutions to the coarse-grained modelling of cascades 4

.1. Model simplification and parameter fitting
The expressions of the cascade output, x Ã n (t), obtained in the previous sections can be used to fit activation data to a small number of parameters. Rather than fitting observations to an entire module of ODEs with n [ N components, the expressions with the gamma function contain three parameters (a (n) , b, n) to describe an optimal cascade, and possibly other real parameters associated with the input (e.g. l for the exponentially decaying input, or v for the periodic stimulus). In particular, note that the first argument of the incomplete gamma function (3.4), which is linked to the cascade length, is a positive real number [19]. Hence when fitting data (figure 2), the estimated length of the cascade is turned into a real-valued parametern [ R, similar to what is done with Hill coefficients to represent multiple mechanistic steps [21].
In figure 2, we present the application of this approach to the fitting of the output of an optimal cascade with two different inputs. We start by generating simulated data from a cascade of length n ¼ 5 with parameters a 1 ¼ 3, . . ,5. One cascade is subject to a constant stimulus R(t) ¼ 1 and the other to an exponentially decaying input R(t) ¼ e Àlt with l ¼ 1. We solve numerically the n-dimensional system of equation (2.2) for both inputs (solid lines in figure 2b,c), and then we generate 'observations' by sampling the output x 5 (t) at times t ¼ f0, 1, . . . , 10g with additive Gaussian noise drawn from N (0, 0:05 2 ). We consider these samples as our 'noisy data' (squares in figure 2b,c) and we fit the gamma function expressions (3.2) 1 and (3.6), respectively, using a Matlab implementation of the squeeze-and-breathe evolutionary Monte Carlo method which is especially appropriate for time course series [20]. 2

Application to near-optimal cascades with random deactivation rates
Strict optimality of cascades [4] requires that all deactivation rates of the proteins be identical (i.e. b i ¼ b for all i). Likewise, our expression for the cascade output in terms of the incomplete gamma function is only strictly valid under the same assumption. Naturally, it is unreasonable to expect identical rates in a biological system. Therefore, we ask the question: if we relax this condition and allow each b i to be an independent and identically distributed (iid) random variable with mean b, can we still approximate the output of the module with a gamma function?
We have tested this idea in figures 3-5. First, we check that cascades with non-identical deactivation rates still achieve maximal amplification when the cascade is of finite length, and we characterize the distribution of cascade lengths observed. Figure 3 shows the histogram of the cascade length at which maximal amplification is achieved for random ensembles of cascades. We consider a step-function input R(t) ¼ 1 with a 1 ¼ 1.2, and we take as a reference an optimal cascade with identical activation rates a i ¼ 1 for i . 1 and deactivation rates b i ¼ b n ¼ (a 1 G) À1=n ¼ 9:6 À1=n , which delivers a gain of G ¼ 8 with an optimal finite length of n ¼ 4 [4]. We then generate 1000 sets of cascades of length n ¼ 1, . . . ,10, with deactivation rates drawn from a distribution b i N ( b n , 0:05 2 ), b n ¼ 9:6 À1=n , i ¼ 1 . . . n and n ¼ 1, . . . ,10 and we record the length at which the maximal amplification occurs. Note that the mean of the deactivation rates depends on the length of the cascade. As shown in figure 3, near-optimal cascades (with normally distributed b i with mean b n ¼ 9:6 À1=n ) achieve maximal amplification for lengths between n ¼ 3 and 5 in 60.4% of cases.
To test whether we can use the gamma function to estimate the parameters of cascades in which the deactivation rates are not identical, we simulated 1000 cascades under a stepfunction input R(t) ¼ 1, with a (n) ¼ 3, n ¼ 5, and random deactivation rates b i N (2, 0:05 2 ). In each cascade, we fitted the parametersâ (n) ,b,n in equation (3.2) to the 'observed' x Ã 5 (t). Figure 4 shows the histograms of the fitted parameters for the 1000 random cascades. The fitted parameters are close to their 'true' values, with the distributions ofâ (n) andn peaked close to their 'true' values, and the distribution of estimated deactivation ratesb normally distributed around its 'true' value.
To check that the outputs for (random) near-optimal cascades can be well approximated using the gamma function expressions, we show in figure 5 that the distribution of asymptotic values of an ensemble of cascades governed by (2.8) with x · * n-1 = a n-1 x* n-2 -b x* n-1 x · * n = a n x* n-1 -b x* n linear n-cascade, up to n + 2 parameters Figure 2. Simplification of a linear activation cascade and fitting with incomplete gamma functions. (a) Schematic of an optimal linear cascade (2.2) and its corresponding equivalent output function (3.2) under a step-function input. The output of the cascade, x n , relays the signal to downstream components of the pathway. Whereas the full n-dimensional model of the cascade has up to n þ 2 parameters (a i , b, n), the condensed expression for the output has three parameters a (n) , b, n. Fitting time-courses of a cascade with two different inputs: (b) a step function and (c) an exponentially decaying stimulus. In both cases, we considered an optimal cascade with n ¼ 5 components and parameters a 1 ¼ 3, a i ¼ 4 for i ¼ 2, . . . ,5, and b ¼ 3. The step-function input was R(t) ¼ 1, t ! 0 and the exponentially decaying input was R(t) ¼ e Àlt with l ¼ 1. The solid lines indicate the solutions to the full system of n ODEs. The squares are 'noisy data' generated from the full model: x 5 (t) sampled at t ¼ f0, 1, . . . , 10g with additive Gaussian noise with standard deviation s ¼ 0.05. The dashed lines are fits of the noisy data using the corresponding incomplete gamma function expressions, equations (3.2) and (3.6). The fits were carried out using the squeeze-and-breathe algorithm [20]. (Online version in colour.) rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160409 b i N ( b n , s 2 ) matches the distribution obtained from the gamma function representation (3.2) with b i N ( b n , s 2 =n). Hence, the gamma function form can be used for near-optimal cascades with random variability in the deactivation parameters, by scaling the variance of the deactivation rates by the length of the cascade (figure 5c).

Cascade reordering: lumped representation of identical blocks
As another deviation from strict optimality, we examine how the output of a weakly activated cascade is modified when a single protein in the cascade has a different deactivation rate. For instance, Chaves et al. [4] considered an auxiliary protein with different deactivation rate at the end of the cascade. We study the effect of such a 'perturbation', and the effect of the position of the perturbation in the cascade. Consider a cascade of n proteins with activation rates a j and deactivation rates b j ¼ b, 8j = i, and b i ¼ b þ 1 for a given node i. First, note that from the Laplace transform of x Ã n (t), it is clear that the position in the cascade of the protein with distinct deactivation b i does not affect the final output This fact allows us to reshuffle the equations of linear cascade models, grouping the blocks with identical deactivation rates, which can thus be lumped upstream in the cascade and replaced with the incomplete gamma function representation. The equations of the perturbed proteins can be placed downstream and take the gamma function of the lumped block as an input. Such reordering can be used to reduce and simplify the model of a cascade without altering the dynamics or timescales (figure 6a).
More explicitly, suppose we have an 1-perturbed cascade of (n þ 1) proteins reordered so that the first n proteins all have deactivation rate b and the (n þ 1)th protein has rate b þ 1. For a step-function input R(t) ¼ 1, t ! 0 we use equation (3.2) to summarize the first n equations, and the equation for the perturbed (n þ 1)th protein becomes then This equation can be solved analytically to give where we have assumed the initial condition x Ã nþ1 (0) ¼ 0. Likewise, when the input is exponentially decaying R(t) ¼ e Àlt , we have that When the initial condition is x Ã nþ1 (0) ¼ 0, the analytical solution for b = l is " # :

ð4:6Þ
We illustrate these points in figure 6. Figure 6b shows the time course of a cascade with six proteins in which the deactivation of the third protein is perturbed. We then reorder the equations so that the perturbed one lies at the bottom. Figure 6c shows the output of the first five reordered equations, given by the gamma function expression (3.6) (dot-dashed line), and the analytical solution of the perturbed protein (which is now the output of the cascade, solid line), given by equation (4.5). Note how the time-courses of the fifth protein in the original and rearranged cascades are different, yet the time course of the sixth protein is identical in both cases, as per our solution. Given the results for random cascades presented above, this approach can be applied to lump sub-cascades of proteins with similar deactivation rates which can then be described compactly through their corresponding gamma function modules.

Simplified modules for activation cascades and delay differential equation models
Experimental observations in signalling cascades are typically concerned with the amplification, distortion and delay introduced in the output. As discussed above, when using ODE models, delays are usually incorporated through the addition of extra equations (and their corresponding extra variables and parameters) corresponding to unmeasured, hidden components, steps or processes in the cascade [22]. This approach can lead to large (high-dimensional) models with many unobservable variables and high numbers of parameters to be identified or fitted [23,24]. Alternatively, modellers often use DDEs to account for the lag between an event and its effect   [25][26][27]. In a DDE, the activity of a variable depends on the state of the system a time t in the past: where the parameter t ! 0 is the delay. Although linear systems of DDEs can in principle be solved analytically using infinite series involving the Lambert function [28,29], such solutions are often impractical to use. We have checked that our results can be applied to model simple delays in linear activation cascades, leading to concise ODE models that capture the delay through the gamma function terms without the need to rely on DDEs (figure 7a). As an example, consider a system with delay modelled with the linear DDE: dp 1 dt ¼â Àbp 1 and dp 2 dt ¼âp 1 (t À t) Àbp 2 : 9 > > = > > ; ð4:7Þ Figure 7b(i) shows the simulated time course ofp 2 (t) (solid line) whenâ ¼ 2,b ¼ 3 and t ¼ 2 with initial conditionŝ p 1 (0) ¼p 2 (0) ¼ 0. This series was numerically obtained with the dde23 solver in Matlab. We then generate our 'observed data' by samplingp 2 at various time points and adding observational random noise from a distribution N (0, 0:05 2 ). We then fit this noisy data to our gamma function expression (3.2): and we estimate the corresponding parameters. Figure 7b(i) shows the fit, as obtained with the squeeze-and-breathe algorithm [20], with estimated parametersâ (n) % 2:270, b % 7:530 andn % 22:107.
To explore the connection between the parameters of the DDE and the best-fit activation cascade model, we simulate the DDE (4.7) with parametersâ ¼ 2 andb ¼ 3 for different values of the delay t [ ½0, 5 and fit models as above. The dependence of the fitted parameters and t is shown in figure 7b(ii). Reassuringly, the amplification factorâ (n) remains relatively constant, whereas the ration=b grows linearly with t. This can be expected from the simple argument that the time delay t in the DDE should be related to the accumulated time needed to traverse n sequential steps with duration 1/b. Hence, a DDE with delay t can be approximated with a linear cascade, by tuning both the length and the deactivation rate of the cascade, i.e. t n=b À 1.  In Beguerisse-Díaz et al. [30], we have used the approach described here to introduce delays in the antioxidant responses of guard cells to abscisic acid and ethylene stimuli during stomatal closure in an ODE model.

Discussion
In this work, the classic model of activation cascades in the weakly activated regime [1] has been re-examined. We have considered the important case where all deactivation rates of the components of the cascade are identical, which was shown to provide optimal amplification in Chaves et al. [4]. Our results show that the output of optimal cascades can be represented exactly by lower incomplete gamma functions, and we show numerically that even when the cascades are near optimal (i.e. when the deactivation rates are iid normal random variables) a gamma function can summarize the cascade by an appropriate rescaling of the parameters. We also show that the position of a protein in the cascade does not affect the final output, so that blocks of proteins with identical deactivation rates can be lumped and represented with incomplete gamma functions. These results allow the reduction of the number of equations and parameters in ODE models without affecting the dynamics or the time scales of the system. We have also shown that in some cases incomplete gamma functions can be used to model delays within systems of ODEs, as an alternative to DDEs. Beyond its application to enzymatic activation cascades, similar mathematical models of cascades could be helpful for the parametrization and modelling of multi-step transcriptional processes, an area of active research in systems and synthetic biology [16,[31][32][33]. In general, model reduction of systems of differential equations remains a challenging and active area of research [34][35][36]. Some methods reduce network models (or modules) based on the topology, effectively finding a minimal kernel that preserves some aspects of the dynamics [37]. Yet, by only considering the topology of the system such methods cannot be guaranteed to preserve time scales or behaviour [38], and are best suited for Boolean models. As Beguerisse-Díaz et al. [30] show, time scales and transients can be crucially linked to the behaviour of a model and cannot be ignored in many cases. Our work introduces a simplified, compact description that can serve to consider delays in ODE models for systems and synthetic biology, and to fit data from experimental observations. Data accessibility. No new data was generated during the course of this research.
Competing interests. We declare we have no competing interests.