HillTau: A fast, compact abstraction for model reduction in biochemical signaling networks.

Signaling networks mediate many aspects of cellular function. The conventional, mechanistically motivated approach to modeling such networks is through mass-action chemistry, which maps directly to biological entities and facilitates experimental tests and predictions. However such models are complex, need many parameters, and are computationally costly. Here we introduce the HillTau form for signaling models. HillTau retains the direct mapping to biological observables, but it uses far fewer parameters, and is 100 to over 1000 times faster than ODE-based methods. In the HillTau formalism, the steady-state concentration of signaling molecules is approximated by the Hill equation, and the dynamics by a time-course tau. We demonstrate its use in implementing several biochemical motifs, including association, inhibition, feedforward and feedback inhibition, bistability, oscillations, and a synaptic switch obeying the BCM rule. The major use-cases for HillTau are system abstraction, model reduction, scaffolds for data-driven optimization, and fast approximations to complex cellular signaling.

This is a PLOS Computational Biology Methods paper.

Introduction

John von Neumann’s elephant haunts mechanistically detailed models. von Neumann was reported to have claimed that he could fit an elephant with 4 parameters, with the implication that models with many parameters are under-constrained and over-fitted [1]. There are two major arguments to exorcise this elephant: that mechanistic detail is needed to address certain kinds of questions; and that in the era of big data it is both easier and less biased to simply build up detailed models with all the available pieces. Here we describe a model formalism, the HillTau form, to help navigate between biological mechanisms and big data on the one hand, and the desirability of condensed model representations that expose the key principles of system function.

Cellular, and particularly synaptic signaling, is notoriously complex. There are an estimated 1400 protein species localized to the postsynaptic density alone [2]. These support a range of functions including synaptic transmission, maintenance, plasticity, activity-driven protein synthesis, metabolic control, and traffic [3].

Mass-action chemistry is a common denominator for mechanistically inspired modeling of these phenomena. This has the key virtue of defining specific biological entities (molecules) and processes (reactions) that map directly to experimental observables. Many studies are based at this level [4,5]. Further levels of mechanistic detail include reaction-diffusion, stochastic chemistry mesoscopic stochastic methods with trapezoidal or cubic meshes [6,7] and even single-particle reaction-diffusion calculations [8,9]. Note that the additional mechanistic detail comes at a considerable computational cost.

A few studies have found ways to lessen the level of detail, typically by focusing on interactions without dynamics (e.g, [10]) or on dynamics with highly reduced interactions (e.g., [11]). Model detail may also be abstracted out through model reduction, which starts from a detailed (usually mass-action or Michaelis-Menten ODE form) model and strips it down to core interactions needed to account for model behavior. Another reduction approach is to identify ‘fast’ reactions in the system, which settle much faster than the overall system, and can be replaced with algebraic relations [12,13]. These are a subset of general approaches to model reduction using quasi-equilibrium and quasi-steady state methods (reviewed in [14]). There are serveral other model reduction techniques (reviewed by Snowden [15]). Most of these methods retain the chemical kinetics formalism using ordinary differential equations (ODEs) to represent mass-action chemistry.

Biochemical signaling models frequently suffer from incomplete parameterization. Thus ‘detailed’ models of signaling pathways, which are of course essential for many kinds of mechanistic analyses and design of experiments, are often under-constrained. In this context, a reduced model is preferable as it requires fewer parameters. One frequently used form specifies rate of change of concentration of each molecule as a weighted sum of input molecule concentrations, which may be passed through a sigmoid to achieve saturation [16-18]. This form is quite similar to neural network models. Thus it lends itself to machine learning approaches to obtain parameters from systematic experimental time-series measurements [17]. The authors obtained relatively sparse interaction weight matrices, thus keeping down the number of parameters. While this formulation is effective at modeling dynamics of molecules in reaction networks, the resultant interaction matrices do not map directly to reaction pathways. Similarly, other formal approaches to model reduction yield very compact models, but the mapping to experimental observables may be quite indirect [15]. Hence it is useful to have a compact chemically-inspired formulation to serve as the core for the model reduction while remaining easy to parameterize and predict using the same quantities that are measured in experiments [19-21]. Indeed, a compact model with few parameters is arguably a better starting point to understand complex signaling with insufficient data, than is a mechanistically detailed model.

Savageau and colleagues have developed the Design Space Toolbox to facilitate a systematic approach to developing reduced signaling and transcriptional network models with specified properties such as multistability [22]. They cast mechanistic models into a Generalized Mass Action form, and this is then analyzed to realize the required phenotypic repertoire. While this is an effective way of obtaining models with desired multi-state properties, it differs in objectives from our goal of having a reduced, very efficient representation of dynamic responses of complex reaction networks such as synaptic signaling.

Efficiency is a specific constraint in developing models of synaptic signaling. On the one hand, many neural functions depend on the nuances of signaling. For example, network properties are quite sensitive to different plasticity rules [23], neuromodulators [24], and mutations [25]. Network models also are expanding to include diffusible messengers controlling cellular activity and blood flow [26]. At the single-cell level, explorations of receptor insertion and clustering [27,28], sequence recognition [29] and synaptic tagging [30,31] all require some level of reference to the chemical signaling. The crux of the problem arises when these studies need to scale beyond one synapse to whole-neuron (up to 104 synapses, [29]) or even network scales (e.g., 109 synapses [32]). Clearly, efficiency in memory and computations is important for such models.

The HillTau form addresses several key concerns with modeling of complex signaling networks. It utilizes only those observable states specified by the user to map directly to the chemistry, thus supporting sparse models that are easier to constrain with limited data. This requires very few parameters, yet behaves similarly to chemical cascades involving multiple intervening steps. Since the user specifies their chosen observables, each can be related directly to observations of concentration over time. The models are small and calculations are highly efficient, being closed-form and event-driven.

Results

We first provide an overview of the HillTau algorithm. Then we illustrate its use to approximate increasingly complex reaction networks. We then show how one can reduce a mass-action model to its HillTau equivalent, with a tradeoff of greater complexity for better accuracy. Finally we carry out some benchmarks of several reduced HillTau models against the original ODE-chemical kinetic models run on two simulators, MOOSE [33] and COPASI [12], and show that HillTau is orders of magnitude faster.

Overview of HillTau algorithm

The name HillTau comes from combining the Hill form for concentration-dependence of a reaction, and tau, the time-course for settling to steady state. In brief, a ‘reaction’ in HillTau uses the Hill equation with modifiers to estimate steady-state values Y∞ of the product of one or several chemical reactions having an input reagent Y, and a Ligand L, with order n:

It may also optionally have a modifier M, with order h:

A modifier changes the effective KA of a reaction, and is controlled by two terms. Kmod determines the half-max concentration of the effect of the modifier. Amod determines what effect the modifier has on the reaction. If Amod< 1, the modifier is inhibitory, else it is excitatory [34]. The steepness of the effect of the modifier is controlled by its order, h.

The HillTau formulation of a reaction also incorporates τ, the time over which the system exponentially approaches this steady-state. We allow for different time-courses τ and τ2 when the concentration is rising or falling:

This exponential form is a good and efficient approximation to the differential equation form for reaction rates (Eq v), so long as the timestep Δt in Eqs iii and iv is smaller than τ (See methods):

The set of elementary HillTau reactions are illustrated in Fig 1, and the details of the calculations are provided in the Methods section.

Fig 1

The HillTau formulation and representation of elementary chemical reactions by a single HillTau reaction.

A: Principle of HillTau formulation. Left: steady-state output values for different levels of the input molecule, computed by the Hill equation. In all simulations in this figure, two input values are used: first 1 μM (red dot) and later 0.2 μM (blue dot). Right: The simulator starts from the current value of the output, and computes the approach to the steady-state as an exponential time-course. Note that these are algebraic calculations, not numerical integration. In this example the output rises from zero toward the red dotted line for 2 seconds. Then the input is changed to 0.2 μM, and now the simulator approaches the steady-state value for this (blue dotted line) with an exponential time-course. B: Key section of the JSON code defining this reaction system. C-G: Inputs (blue) and simulated time-course of outputs (orange) for seven different reactions. Each is represented by a similar HillTau reaction but with different parameters (see S1 Data). In all cases the HillTau output onset is identical to the output computed using numerical integration of a single reaction expressed as mass-action chemical kinetics. Decay time-course may differ from onset time-course in mass-action. C: Binding. D: 2nd order binding. E: Conversion. F: Inhibition, conceptually equivalent to removal of output molecules by binding of input to the output molecule, and sequestration of the resultant complex. G: Variant of binding reaction, in which there is a fixed baseline of 0.5 μM, and the system has different on (tau = 1s) and off (tau2 = 5s) time courses. H. Same as reaction 1, but with a modifier term that strengthens input affinity. I. Same as reaction 1, but with a gain term that multiplies the output, in this case by a factor of 2.

The motivation for this formalism is that the steady-state value of a cascade of binding reactions, or of enzyme reactions with a fixed rate back-reaction, can be approximated by a Hill function (Methods). Further, the time-course of approach to steady-state is typically governed by the slowest reaction, and this can be approximated as an exponential settling function (Methods).

Note that we do not assume that the input, activator and modifier act in a single mass-action chemical step. Indeed, HillTau is most effective for model reduction when one can fit several mass-action steps using one HillTau ‘reaction’.

Since Eqs i to iv are analytic, one can do this calculation in an event-driven manner. HillTau achieves sparseness and simplicity by approximating many steps with a single ‘reaction’, considering only those intermediates that are needed for readouts or for improved precision. It achieves speed because the models are smaller, and by using event-driven calculations rather than numerical integration.

Most reaction networks cascade through many layers of reactions. HillTau evaluates each upstream layer before downstream ones. It first builds a dependency graph of all reactions. This is done by identifying input molecules as layer 0, and successively ranking all reactions that depend only on layer 0 as layer 1, reactions that depend on layers 0 to 1 as layer 2 and so on.

HillTau identifies feedback loops by reactions which do not resolve into the above layers. Based on ordering of reactions in the model definition, it picks a reaction to ‘break’ the loop, and assigns it to layer N+1, where N was the previously deepest layer. It then repeats the process of layer assignment, including further loop-breaking if needed.

During evaluation of a single step in HillTau, all the steady-state and time-course calculations are completed for layer 1, then layer 2 is calculated, and so on. Thus each layer receives the inputs appropriate to the current time before doing its evaluation. In cases where update events are separated by periods greater than the shortest τ in the system, additional time-steps are inserted to maintain accuracy (Methods). For typical use-cases, such as synaptic plasticity models, the event interval is shorter than the time-courses in the model (typically ~1 sec) and hence only a single step is taken. In cases where HillTau inserts additional time-steps for accuracy, it is done behind the scenes of the same event-driven programming interface. If there is feedback, then again one has to use event intervals shorter than the shortest τ in the feedback loop. A factor of 10 usually gives good convergence (e.g., S7A and S7D Fig).

In summary, HillTau uses analytic evaluation of reaction outputs based on a Hill-like form and exponential settling, and propagates the evaluation through successive layers of the reaction network for each event time. Events can be stimuli or points in a time-series for sampling system time-evolution.

The HillTau form can model a range of chemical signaling motifs

We implemented a range of elementary chemical signaling functions to illustrate the use of the HillTau form (Methods, Fig 1). The HillTau versions of most of these reactions have an exact fit to their mass-action counterparts (S1 Fig). We further implemented key signaling motifs, including feedback inhibition, oscillation, and bistables (Fig 2). To do this, we constructed minimal HillTau schemes that incorporated the essential elements of each of these motifs. We developed an optimizer program mash.py (MASH: Model Abstraction from SBML to HillTau, see Methods) to tune parameters of the HillTau models to match the outputs to the original mass-action or ODE versions. MASH runs the reference model through a range of stimuli designed to explore its input-output properties, and then uses numerical optimization methods from scipy.optimize to tune parameters so that the HillTau model produces a good fit to the original. We used normalized RMS difference between the traces as a measure of goodness of fit. In Fig 2A–2C we compare feedback inhibition implemented in mass-action (5 reactions, 7 species, 2A), HillTau (2 reactions, 3 species, 2B), and run for a square pulse input (2C). The feedback inhibition model is well approximated by the HillTau version to within 4% normalized RMS deviation.

Fig 2

HillTau models of key signaling motifs.

A-C: Feedback inhibition. A: Mass-action reaction scheme for feedback inhibition, involving 7 molecules and 5 reactions. B: HillTau version. Each box represents a molecule. If there are input arrows to the box it means there is a reaction whose product is the named molecule. Input arrows can be either inputs (reagents), activators, inhibitors, or modifiers. This reaction consists of 3 molecules (input, fb, and output) and 2 reactions (fb and output). C: Simulations for mass-action (blue) andHillTau (orange) versions of feedback inhibition. The green trace is the input molecule. D-F: Oscillator from ultrasensitive MAPK cascade, taken from [35]. D: Output of simulation. Blue is ODE output and orange is HillTau. E: ODE model. This uses 15 molecules, and 11 reactions. MAPK-pp is the molecular species used as output of the oscillator. F: HillTau reaction scheme for oscillator, using 5 molecules and 3 reactions. The concentration of the ‘output’ molecule is plotted. G: HillTau model of bistable system, involving 4 molecules and 2 reactions. H: Phase plot showing stable states of system as the intersection points between the steady-state dose-response curves. This was generated by varying the feedback molecule fb, and measuring output (brown curve), and then varying the output molecule and measuring fb (pink curve). I: Time-series illustration of state switching in the bistable. As before, output is in brown and fb in pink. The Y axes of H and I are the same to show that the steady-state output levels (brown) match. The system starts in the low state. At 20 s a small excitatory input stim is given which fails to switch the state. At 40 s a strong input causes switching to the high state. At 60 s a weak inhibitory input fails to turn it off, but at 80 s a strong inhibitory input returns the state to baseline. Excitatory and inhibitory inputs were delivered by transiently setting the level of stim to high or low values.

Next, we implemented a HillTau version of a mitogen-activated protein kinase (MAPK) feedback oscillation model having 11 reactions and 15 species, Fig 2E [35]. We used three HillTau reactions to map to the key components of the original ODE model. First, we used a reaction to represent the basic MAPK cascade. Second, we provided an output reaction to represent the phosphorylation of the MAPK molecule by the cascade. While it was possible to use this output signal to inhibit the cascade, we found we had to implement a separate reaction for the negative feedback step to introduce a longer delay to match the observed oscillations.

Having constructed the model structure, we next fit the HillTau model to the original ODE model using MASH. As initial parameter estimates, we used taus of the order of the oscillatory period, and KA of the same order as the (known) molecular concentrations. We first fit the initial output transients. Then we ran it for a complete cycle. Finally we stretched the fit time to include a few cycles. This incremental increase in time was necessary because our RMS scoring function gives very poor scores for otherwise good models if a phase mismatch builds up over a few cycles. The final reduced HillTau version (3 reactions, 4 species, Fig 2F) had similar period and amplitude (Fig 2D), and it fit the waveform to within ~7.6% normalized RMS deviation.

Finally, we made a HillTau version of a chemical bistable switch using just 2 reactions (Fig 2G). We demonstrated that the HillTau form works with the standard dose-response (null-cline) approach to estimating steady states (Fig 2H) and showed that the resulting switch exhibits high and low states that are triggered by transient inputs (Fig 2I).

Thus the HillTau form can efficiently represent a range of important signaling motifs and their dynamics, including feedback inhibition. oscillations, and bistability.

The HillTau form compactly represents bidirectional synaptic plasticity

Synaptic plasticity is one of the most-modeled neuronal signaling processes [4,36]. The key features that have been represented include stimulus strength-dependence, timing dependence, and long-term state storage [37]. A few studies have come up with rather detailed models to implement each of these processes [31,38,39]. As an illustration of all these properties in the HillTau system, we implemented bidirectional synaptic plasticity including long-term synaptic state changes (Fig 3). One of the interesting aspects of synaptic plasticity is that in many systems, the same input modality (typically read out as Ca2+ concentration) can give rise to both synaptic depression and potentiation. This has significant theoretical implications and an abstract rule for this bidirectional plasticity was proposed by Bienenstock et al. (the BCM rule, [40]). We first devised a simplified mass-action version of the BCM rule using 9 molecules and 6 reactions (Fig 3A). The species p_AMPAR is the phosphorylated form of the receptor, assumed to be inserted into the synapse. Here, resting Ca2+ does not alter the state of the model; low Ca2+ causes depotentiation (that is, reduction of receptor levels), and high Ca2+ causes potentiation (Fig 3C, 3D and 3E). We then implemented a BCM model in just 3 reactions in HillTau (Fig 3B). We used the program mash.py to fit the HillTau model to the reference mass-action version using a set of generic time-series and dose-response stimuli (Methods). We obtained a normalized RMS fit of ~2.3%. When we used the fitted model for Fig 3, we obtained fits of ~5.2%, 8.9% and 2.3% for panels C, D and E respectively even though the model had not been tuned to these stimuli. As a further elaboration, we introduced a bistable switch for long-term retention of synapse state, which was driven bidirectionally by the BCM rule (Fig 3F). The bistable switch, derived from Calcium-calmodulin Type II kinase (CaMKII) signaling, controls receptor insertion. Using this model we delivered a typical potentiating stimulus (strong but brief Ca2+ input), leading to sustained synaptic AMPAR elevation. We followed this with a typical long-term depression stimulus (modest but sustained Ca2+ input), which turned the switch off again and led to reduction in AMPAR (Fig 3G). This composite model required 4 reactions and one summation function in the HillTau form. Several mass-action models of synaptic state switches include these elements (e.g., [4,36,38,41]) and they typically involve far more molecules and reactions (e.g., the Hayer and Bhalla 2005 model used 133 molecules and 215 reactions [38]).

Fig 3

HillTau version of synaptic plasticity rules.

A. Mass-action model for generating Bienenstock-Cooper-Munro (BCM) rule for synaptic plasticity. p_AMPAR is the phospho-receptor, and is the output of the model. It is assumed to localize to the synapse and is thus also referred to as synAMPAR. Calcium triggers both an inhibitor (Calcineurin, CaN) and a stimulus (CaMKII) for receptor phosphorylation and insertion into the synapse. CaN activates at lower [Ca2+], so there is initially a reduction in p_AMPAR. CaMKII is present at very high levels, so at higher Ca2+ it out-competes CaN to give an increase in p_AMPAR. B. BCM rule implemented in HillTau. Here the species synAMPAR is the output of the model. C-E: Comparison of mass-action model p_AMPAR with HillTau model synAMPAR. Orange is HillTau, blue is mass action. C: 1s stimulus at 0.5 μM Ca2+ gives a reduction in synaptically localized AMPAR (synAMPAR). D: 1s stimulus at 5 μM Ca2+ gives an increase in synAMPAR. E: Dose-response curve of steady-state synAMPAR as a function of [Ca2+] for mass-action (blue) and HillTau (orange) models. In both cases settling time for each point was 1000s. F: Schematic for BCM rule model feeding into bistable model, implemented in HillTau. The circular node labeled Σ represents weighted summation of multiple inputs. G: time-course of simulation of bidirectional plasticity using different Ca2+ stimuli. At t = 20, a 1s stimulus of 2μM Ca2+ (green trace) causes a transition to the active state, using synaptic AMPAR (maroon trace) as a readout. At t = 50s, a 30s stimulus of 0.3 μM Ca2+ pulls the system back to resting state.

Overall, these examples illustrate how compact HillTaumodels can represent both the bidirectional induction of plasticity, and also long-term maintenance of synaptic state.

HillTau models can be optimized to fit biochemical measurements

The above examples illustrate how HillTau can represent biological signaling motifs, and build them up into networks with interesting computational properties. We next approached a complementary problem in signaling, namely, to take a complex signaling system, and fit simple HillTau models to it. This provides a way to perform model reduction and to infer computational properties. The basic flowchart is illustrated in Fig 4. This flowchart addresses both the heuristics of defining model topology, and of parameter fitting.

Fig 4

Flowchart for model building using HillTau.

Left: flowchart. Right top: initial reaction with inputs A, B, and C. Right middle: successive local increments to model, introducing reactions X and Y respectively. Right lower: Final HillTau scheme with good fit at all stages.

The heuristics for defining model topology are as follows.

Identify inputs and key readout molecules. These readouts may be important (and experimentally measured) intermediate signaling molecules in a reaction network, or the end-products of a cascade.

Assign a reaction for each readout molecule, with an input as an upstream substrate or inactive state of the molecule, an activator (or inhibitor, see methods) and optionally, a modifier. Together these control the level of the readout molecule.

In case a molecule has multiple inputs, bring in additional reaction steps based on the known reaction mechanisms. For example, if we have BDNF, EGF and Ca all controlling ERKII activity, then we could specify an intermediate step where the two receptor tyrosine kinase ligands converge, and this combination is an activator for the ERKII reaction with Ca as a modifier.

In case a readout is simply the sum of multiple active states of a molecule, use an equation to define this summation.

Obtain best model fit as per flowchart in Fig 4. If model accuracy does not meet criteria for your objectives, identify poorly performing intermediate readouts and insert further intermediate reaction steps.

Note that in principle many of these steps can be automated. For example, one can generate a family of models algorithmically (e.g., [42]) and optimize over topology as well as parameters, but this is out of the scope of the current study. Other algorithmic approaches for model reduction are discussed in [14].

In the current paper we have used the above heuristics to generate HillTau schemes (model topology) by hand.

There are also simple steps to obtain initial parameter estimates for each HillTau ‘reaction’:

KA from the activator concentration at half-maximum of the experimental activation curve, or directly from mass-action model rates.

Time-course τ from the experimental time-course of the reaction, or from the slowest intermediate step of a mass-action model. If there is a distinct time-course when the reaction turns off, use this as τ2.

When the modulator is present, assign Kmod and Amod from the half-maximum and the steepness of modulation curve.

This approach works in the same way for model construction from experimental response curves, and for model reduction from response curves taken from detailed models. Following generation of an initial, roughly parameterized HillTau model, we can deploy the model fitting approaches described in FindSim[43]. In brief, FindSim provides a Python-based framework for matching models to experiments. It codifies the experiment design (e.g., time-series, dose-response, bar-chart) and experimental results into a single machine-readable file. FindSim runs the experiment on the model and returns a numerical score for goodness of fit. The model may be defined in SBML (run using MOOSE) or using HillTau. Thus FindSim can be used as the scoring function for optimizing model fit to experiments using a variety of optimization methods available in scipy.optimize.

For the special case of model reduction, where we already have a detailed SBML model and wish to fit a reduced, HillTau version, the utility MASH provides a shortcut alternative to the FindSim and optimization pipeline (discussed above and in the Methods section).

As an example of this flowchart and the use of HillTau fitting to match an existing, detailed chemical ODE model, we derived a HillTau model of synaptic activity-triggered protein synthesis. Our reference data was obtained by running a series of ‘experiments’ on a published model implemented in mass-action kinetics [44]. The original model was based on numerous experiments, and included 123 molecules and 120 reactions (Fig 5A). The input pathways were Ca2+ and brain-derived neurotrophic factor (BDNF), and the final output was protein synthesis rate.

Fig 5

Model fitting and model reduction.

A: Block diagram of source model with 123 molecules and 120 reactions, from [44]. B: First pass reduced model in HillTau, with 1 reaction and 4 molecules. C: Final reduced HillTau model including activated S6K (aS6K) and activated CaMKIII (aCaMKIII) as intermediate readouts which also were fit to data. This model has 4 reactions and 10 molecules. D-K: Eight ‘experiments’ on reference and HillTau models, not part of stimulus set used to tune parameters for HillTau version. In all cases protein production rate is readout. Blue plots are reference, orange areHillTau. D: BDNF@3.7 nM + Ca2+@0.2 μM, 900 seconds. E: BDNF@3.7nM, Ca2+@1μM. F: 3 pulses of BDNF@3.7 nM for 5s, coincident with Ca2+@10μM for 1s, pulses separated by 300 s. G: Same as F, but Ca2+ held at baseline of 0.08 μM. H: Dose-response of protein vs. Ca2+, holding BDNF fixed at 3.7 nM. I: Dose-response of protein vs BDNF, holding Ca2+ fixed at 0.08 μM. J: Protein production rate for fixed BDNF@3.7 nM, where Ca2+ was given in 1 second pulses at intervals of 300, 120, 60 and 10 seconds; each pulse sequence lasting for 1200 s.. K: Dose response of protein production rate vs. Ca2+, holding BDNF at basal levels of 0.05 nM. The average normalized RMS difference across the eight ‘experiments’ was under 10%, and in all cases the qualitative properties such as direction of change, matched well.

We started with the most reduced form, a single reaction to replace the entire synaptic protein synthesis network. We specified amino acids as the input, BDNF as an activator, Ca as a modifier, and protein as the product of this reaction, (Fig 5B) to obtain our starting reduced model. We used MASH to carry out the optimization (Methods). In a model with a single reaction, MASH obtained a fit of about 11% normalized RMS. This is remarkable for such a simple model. It does well with the dose-response experiments (Scores of 7% and 5%, S2 Fig). However, it does not do a good job of replicating the dynamics of the experimental data, achieving scores in the range of 10% to 42%. (S2 Fig). We therefore increased the HillTau model detail. To do this, we introduced two additional key intermediates into the HillTau model: S6K, and CaMKIII. We first made a HillTau model involving inputs to S6K alone. Based on the known pathway, we chose BDNF as the activator, and Ca as the modulator for S6K. The S6K responses to combinations of these two inputs were used in MASH to obtain a fit to within 3.3% (Individual panel fits were between 3% and 22%, see S3 Fig). We then held S6K parameters fixed while we fit CaMKIII. CaMKIII is activated by Ca, and modulated by S6K. MASH gave a CaMKIII fit of 1.9% (Individual panels 2% to 20% but most of the poor scores were small differences at baseline; S4 Fig).

Collectively, the optimizations for S6K and CaMKIII correspond to the inner loop of Fig 4. Finally, we added a final protein synthesis reaction that took the already fitted S6K and CaMKIII activity as activator and modifier. We held the earlier reactions (S6K and CaMKIII) fixed, and optimized only the protein synthesis reaction. After this, the composite model fit the optimization waveform in MASH to within 3.2% (normalized RMS), and the figure panels fit within a mean of 9.3% (Fig 5D–5K).

At this stage one could choose to perform further optimization in a couple of ways. We could have obtained closer fits had we optimized to the same stimuli as in the figure panels. Instead we used more general input time-series and dose-response stimuli to the MASH optimizer to see how well the model would generalize. This gives a HillTau model that behaves well across a wider range of conditions than the experiments in Fig 5. Second, there is a small systematic difference in baseline in panels 5F and 5G arising from a difference in output at resting BDNF, seen in panel 5I. The introduction of additional intermediate reactions in the model as per Fig 4 could further improve the fit. Such tradeoffs between generality, accuracy, and model complexity are common, and given that HillTau is meant for building compact models we considered this fit sufficiently good for illustration.

Overall, the abstract HillTau model captures many of the key properties of the mass-action system. These include steady-state and time-series responses of two inputs (BDNF and Ca), and three readouts: S6K, CaMKIII, and the end-product protein (Figs S3, S4, and S5 respectively). This is a highly effective dimension reduction, from over 360 to 31 parameters.

It is important to note that the efficiency of HillTau made the optimization calculations quite tractable. In the final optimization run for the entire model, the single reference ODE run took ~60 seconds, and the cumulative time for 746 HillTau evaluations was around 18 seconds. The optimization algorithm itself took about 75 seconds, excluding function evaluations.

Can we create a reverse mapping from these simplified HillTau models to ODE forms? A close but not exact mapping is obtained by taking the small-time limit of the HillTau event-driven form (Equations 3.x) and converting to an ODE (rate) form (Equations 4.x). ODE equations are supported by many systems biology simulators. It is not an exact mapping because HillTau may use different values for rising and falling time-courses (tau and tau2), whereas the ODE form can accommodate just a single value, tau. We implemented this conversion in a program, ht2sbml.py, which is provided on the GitHub repository for HillTau. Using this mapping we were able to export HillTau to SBML, and tested that SBML-capable simulators such as COPASI could run the reduced, ODE form models, and give approximately matching results (Methods, S5 Fig). Thus we can use the HillTau toolchain to make reasonably reduced ODE models, though these are neither as efficient as their HillTau counterparts, nor do they have the same capabilities to use two time-courses to improve model fitting.

In summary, we developed a systematic procedure for developing reduced HillTau models to fit mass-action simulations, including a model optimization utility MASH. We illustrate this procedure by developing a HillTau model of 10 molecules and 4 reactions to fit a mass-action model having 123 molecules and 120 reactions. The resultant HillTau models generalize well and the fit improves when intermediate reaction steps are added.

HillTau models are compact and efficient

We next took a set of HillTau models of various levels of complexity, and compared various measures of computational cost with the ODE equivalents (Table 1 and Fig 6.)

Table 1

Parameters used to define a HillTau model.

Concentration units can be any of M, mM, uM, and nM, and are specified in the JSON file. The default concentration units are mM.

Parameter	Units	Meaning	Default	Required?
ConcInit	Concentration (can be any of M, mM, uM, nM)	Initial concentration. Only applies to species definitions	0	No
KA	Concentration	Association constant	N/A	Yes
τ	Time (seconds)	Time course for relaxation to steady state	N/A	Yes
τ2	Time (seconds)	Time course of relaxation if output is falling.	τ	No
Gain	None	Scaling factor for reaction output. Used to indicate enzymatic amplification.	1	No
Baseline	Concentration	Baseline value of reaction output	0	No
Kmod	Concentration	Half-saturation concentration for modifier	N/A	Only if modifier molecule is specified.
Amod	None	Activation term for modifier	4	No
Nmod	None	Order of modifier action	1	No

Fig 6

Efficiency of HillTau method.

A. Scaling of number of parameters. The HillTau form uses far fewer parameters than mass-action, and becomes relatively more concise with larger models. B, C: run-time scaling. All run-times are expressed in μs of wall-clock time to execute 1 second of simulation time. Calculations were done on an Intel(R) i7-7700HQ processor, running Ubuntu 18.04. B. Comparison of run-times of the HillTau and mass-action forms, where the numerically intensive sections of HillTau were implemented in C++, and C++ was also used for ODE calculations in two simulators: MOOSE and COPASI. Due to the combination of model reduction and efficient calculation, HillTau has a huge efficiency advantage which grows to over 3 orders of magnitude for larger models. The accuracy with which this set of HillTau models fit their ODE counterparts was in the range 3 to 9%. C. HillTau model run-time increases with number of parameters. R2 = 0.79, slope = 0.015, intercept = 0.006, units in μs/s.

We first compared model complexity, measured as the number of parameters needed to specify the model. The number of parameters scales roughly as

# of molecular species + 2 * # of reactions.

This is a slight overestimate, since some of the molecules are state variables and we do not need initial concentration values for them. Here we consider state variables to be those which are computed, as opposed to defined using initial conditions. In ODE models we estimate this by counting the number of rate terms plus the number of molecular species with a non-zero initial value. In HillTau models we count the rate terms and the species which are not reaction outputs. This yields the approximate scaling terms below. Each reaction needs two parameters, Kf and Kb for conversion reactions, and Km and kcat for enzymes. We sampled from among the mass-action models presented in the above sections, ranging from 3 to over 360 parameters, and included an additional model with almost 750 parameters. (Fig 6A). The HillTau form had a similar scaling with molecules and reactions, except that HillTau also allows for a number of optional terms in the reactions so the average scaling is somewhat larger than 2 * # of reactions. We found that the HillTau form became increasingly effective at model reduction for larger models. Note that here the optimization goal was to obtain a single end-point response (3 end points in the case of the model in Fig 5). Further reactions would be needed to also represent intermediate pathway readouts.

We then examined run-time efficiency. We took run-times for two ODE simulators, MOOSE and COPASI, whose numerical cores are in C++. We compared these with matching HillTau models run using the C++ version of HillTau (Methods, Fig 6B). For small models, HILLTAU was typically about 100 times faster than the ODE calculations, but for large models HillTau became over 3 to 4 orders of magnitude faster (Fig 6B). This set of HillTau models fit their ODE counterparts within 3 to 9% accuracy, but note that the approximations were inherent in the HillTau model structure, and did not arise from lack of numerical convergence. The run-time for HillTau models grew with the number of parameters (Fig 6C, slope = 0.015 μs/s, R2 = 0.79), but also had a dependence on model stiffness due to the requirement that the internal timestep should be smaller than the smallest τ in the model. This suggests that the HillTau calculation cost could be further improved by utilization of a variable-timestep similar to methods for ODE solutions frequently used for chemical kinetic calculations.

Dose-response experiments are particularly efficient to compute using HillTau. An inefficient way to do these for ODE models is to run them out to steady-state for each successive dose. This may take a long while especially if the system is stiff or converges slowly. It is also possible to use linear algebraic root-finding to find the steady-state value in one step, possibly following a short presimulation to bring the system closer to the steady-state [45,46]. In HillTau, the form itself incorporates the steady-state value, so in principle one could leap to the steady value in one step. To be more conservative, the HillTau does so in 10 steps to smooth out transients and to allow any feedback signals to propagate through the system. As an illustration, COPASI performed the steady-state calculation for a large model (accession 92, DOQCS, ~100 reactions, [44]) in ~1 second after some relaxation of convergence criteria. The HillTau equivalent model (same as final model in Fig 5) took about 4 microseconds.

Overall, HillTau models are compact and highly efficient compared to ODE-solved mass-action models. The efficiency improves for larger models.

Discussion

We have designed HillTau, a compact, computationally efficient abstraction of chemical signaling that is particularly effective in building reduced models of complex signaling networks. It uses an event-driven algebraic representation based on the Hill equation and exponential relaxation to steady state. HillTau is effective in representing a range of chemical signaling motifs and complex synaptic models, using biological observables of molecules, reactions, association constants and time-courses. We show its applicability for model reduction by optimizing the fit of its responses to those of a reference mass-action model. This generates very compact models. A similar optimization approach works to build a HillTau model directly from experimental data. ThusHillTau addresses many of the concerns of model-building with limited data, and serves as a scaffold for eventual development of more detailed models.

HillTau is phenomenological and semi-heuristic, in that it uses the Hill equation to achieve concentration dependencies that fit well, but ignores many intervening chemical steps. This combination gives it the strong points indicated above, namely speed, compactness, and consistent mapping to experimental observables, but it also sets out clear limitations. Foremost among these is that it can only make limited predictions on detailed pathway chemistry, since it is missing many reaction steps. For example, a HillTau model would be limited in its ability to predict drug targets or side-effects because it may have lumped together potential molecular targets into a single reaction step. It is, however, quite effective in representing and predicting emergent signaling properties because it captures dynamics and topology of signaling networks.

The current HillTau formulation is limited in its handling of two important aspects of signaling in neurons: stochasticity and diffusion. These phenomena are out of the scope of our current implementation, which has focused on development, validation, simplicity and speed. Many biochemical signaling processes experience substantial stochasticity, particularly in small-volume systems such as the synapse which is a target of our modeling. One possible way to introduce stochasticity would be through the linear noise approximation of the chemical Langevin equation [47], which if used in an event-driven manner could be quite efficient. We anticipate it will take extensive validation to establish its utility in the HillTau framework. Similarly, there are potential ways to elaborate upon HillTau to use an event-driven approximation to diffusion, but these will require later follow-up.

Based on these attributes, we discuss four major use-cases for HillTau: model reduction, system abstraction, scaffolds for data-driven optimization, and efficient approximations to complex cellular signaling.

Model reduction

Several algorithmic approaches have been brought to bear on the model reduction problem, including collapsing multiple mass-action steps into one [48], and power-law generalizations of mass-action signaling [16]. Simulators such as COPASI provide utilities for partitioning a reaction system into fast and slow reactions, thus allowing one to approximate the fast steps with their algebraic counterparts [12,13]. In our test model (Accession 92, DOQCS, [44]) we were able to partition about 80% of the steps into the fast domain (<100 s) using the ILDM method implemented in COPASI. With HillTau one can use a well-known heuristic/optimization approach to simplifying large networks [49,50]. This reduced the same test model down to 4 reactions, about 25-fold (Figs 4 and 5). The approach reported in our current study relies on the modeler starting from a minimal input-> output mapping, then iteratively picking relevant major nodes in the network, and optimizing each subset of the model to fit the data (Fig 4).

Thus one can converge on the minimal set of intermediate nodes (illustrated in Figs 4, 5, and S2, S3, S4) to achieve the desired accuracy of model fit to data. Like other model reduction approaches, this minimal set of nodes is a compromise between available data and model accuracy [15]. Unlike several other reduction approaches, HillTau retains a direct mapping to observable biological entities. Indeed, the HillTau representation of a signaling node may be closer to the conventional intuition based on pathway schematics, than is a full mass-action reaction representation. Like pathway diagrams, each HillTau reaction receives excitatory, inhibitory and modulatory inputs. A further point of similarity is the HillTau models may condense several intermediate steps into a single node on the reaction network. A more subtle point of similarity is that pathway block diagrams typically assume implicit back reactions and decay of activity when stimuli are removed. This too is built into how HillTau reactions work. In contrast to these simple mappings from pathway diagrams to the HillTau form, it is often difficult to map between signaling diagrams and the full mass-action reaction schemes [51,52]. While previous model reduction studies have worked on different pathways than the synaptically biased set explored in our study, a comparison with the survey of methods in [15], suggests that HillTau achieves as good or better model reduction for large models than most other methods.

System abstraction

Next, system abstraction and functional modules help to make sense of complex biological signaling. We propose that HillTau forms a useful tool for arriving at functional modules in complex signaling networks. Such modules have long been considered a conceptual basis for understanding complex signaling [53]. Typically they have been ascertained by manual inspection and dynamical analysis of components of signaling networks, for example, the nested feedback loops in the cell cycle [54]. A more scalable approach to uncovering such modules is to use graph theory for motif analysis on detailed mass-action models, but this approach loses key aspects of system dynamics [55]. With the HillTau formalism and our procedures for model reduction, we are able to generate highly reduced reaction graphs that nevertheless support rather accurate dynamics. The formalism encourages models that can be readily mapped to biology. ThusHillTau supports a data-driven approach to arrive at functional modules.

While functional modules are good for analysis, we note that biology does not necessarily partition signaling networks into neat modules [52]. Indeed, cross-talk between pathways is common. HillTau supports explicit cross-talk interactions, but does not introduce implicit interactions. In this regard it differs from mass-action reaction systems, in which mechanisms such as back-reactions or enzyme sequestration may introduce subtle effects, such as implicit feedback. For example, one can achieve bistability through multistage phosphorylation/dephosphorylation cascades [56]. To represent such effects in HillTau one would have to introduce explicit feedback steps between reactions, such as in the bistability example in Fig 2. Similarly, interesting behavior emerging from chemical saturation, such as zero-order ultrasensitivity [57] would require the use of the explicit math expressions supported by HillTau.

Protein-protein interaction networks are a commonly derived form of abstract networks. These can be purely topological, or may also include reaction dynamics [17]. Can these be parsed into HillTau networks? To first order, protein interaction networks lack the distinction that HillTau ‘reactions’ make between inputs, activators and modulators. Additional information is needed to disambiguate these. Data from sources such as pathway maps, protein domain properties, or Gene Ontology relationships are required to resolve HillTau topologies for a given protein-protein network [10].

Next, can the rates be assigned? In networks that include dynamics (e.g., [17]) this is relatively straightforward to accomplish. In Figs 4 and 5 we illustrate how experimental data, or simulated dynamics of an existing model can be used to parameterize a HillTau model. The same approach could utilize the original timing constraints that went into the Nyman model [17]. Alternatively, a program similar to MASH could explore dynamics of the reference (protein-protein network) model, and use the output to search for parameters of the corresponding HillTau model.

ThusHillTau promotes abstraction through model reduction. The abstracted models expose all interactions explicitly.

Scaffolds for model fitting

Third, HillTau is a useful intermediate step, or scaffold, for model fitting of large mass-action models. Direct model-fitting is difficult in at least two ways: there are typically far fewer experiments than parameters, and it is computationally costly to run a large ODE model many times for carrying out an optimization approach to model fitting. We propose that the HillTau form may provide a useful bridge on both these counts. As we have illustrated in Figs 2, 3 and 5, HillTau models lend themselves to fitting to experiments because they have few parameters and they run quickly. Several advantages accrue from an initial pass to make and fit a HillTau model. 1. In building and optimizing a HillTau model, the optimization dataset will be use-tested, and gaps identified. 2. The essential pathway structure of the model will be defined by the HillTau model, and key interactions identified. The mass-action model must, at minimum, incorporate these interactions. 3. The parameters of the HillTau model set bounds for those of the detailed reaction sets. For example, the time-course of any individual mass-action reaction step must be faster than the HillTau reaction in which it is embedded.

Efficient approximations to complex signaling

Finally, HillTau models are useful because they are efficient. One of the key use-cases envisioned for HillTau was to model complex cellular signaling, with synaptic signaling as an exemplar. Several of the examples in the current paper are in this domain, specifically the BCM curve (Fig 3A–3E, the coupled BCM curve with bistables (Fig 3F and 3G), and a synaptic protein synthesis pathway (Fig 5). While even these large ODE signaling models run somewhat faster than wall-clock time (Fig 6), there are at least two cases where far greater efficiency is desirable. First, as mentioned above, model parameter optimization requires a large number of evaluations of complete simulations (100s to 1000s in our experience). To perform a single evaluation, these synaptic simulations may have to run for many thousands of seconds of simulation time to compare with typical plasticity experiments [58]. Further, one typically optimizes a pathway to fit numerous experiments, all of which must be simulated for each evaluation. Together, this is computationally expensive. In our HillTau optimizations using MASH, the total simulation time for large numbers of pathway simulations was typically even smaller than the time spent by the minimizer code itself.

A second use case for highly efficient signaling models is in synaptic signaling. A single neuron may have over 10,000 spines, and there may be many such neurons in a network. If each synapse is to implement a complex biochemical pathway the computational costs may be formidable. Network plasticity models [59], and cellular sequence selectivity models [29] are examples of this scale of model. Indeed, Higgins et al. [59] have used an efficient event-driven calculation of synaptic weights with a similar exponential decay calculation as in HillTau. HillTau signaling provides a way to implement biologically detailed synaptic dynamics in every synapse, even in large networks.

In summary, the HillTau form and its supporting toolkits for running and optimizing models provide a compact, efficient way to perform data-driven abstraction of complex signaling models.

Methods

HillTau formulation

The HillTau formulation is an event-driven variant of a Hill equation with modifiers. It is specified and evaluated in two stages: the steady-state value, and an exponential time-course of approach to steady-state. As detailed below, reactions in a HillTau model are evaluated in successive layers such that the next estimate for steady-state value of a given layer depends only on boundary conditions and on outputs from preceding layers. In the equations below we expand out the equations for a range of use-cases. In equations 1.x we specify the steady-state values for each use case. In equations 2.x we define the approach to steady state. In equations 3.x we combine Eqs 1 and 2 to summarize the evaluations done in HillTau. In equations 4.x we provide interpretations of the HillTau equations as rate terms which can be evaluated by regular ODE solvers and form the basis for the SBML export of HillTau model systems. However, the definitive form of HillTau is event-driven and the rate-term form should be regarded as a convenient but approximate mapping to conventional mass-action solvers. In equations 5.x we provide a motivation for the form of HillTau as an approximation to complex signaling chemistry.

HillTausteady-state

The HillTau formulation stipulates that the steady-state level Y of each signaling step (which may involve multiple chemical steps) is approximated by a Hill function

Here Y is the input concentration to this signaling step, where Yinput is either a molecule with a predefined concentration (i.e, a boundary condition) or coming from an upstream reaction. Likewise, the reactant L can either be predefined or come from an upstream reaction. KA is the association constant of L with Y.

We elaborate this slightly to accommodate an optional gain term:

We utilize a slightly modified form to permit the Ligand L to act in an inhibitory manner:

In cases where there are modifiers, we include a further term based on the analysis of Hofmeyer and Cornish-Bowden [34]:

Here M is the concentration of the modifier, Kmod is the half-effect value, Amod is the modifier action, and h is the order of the modifier. The modifier acts in an inhibitory manner when Amod< 1, and as an activator when Amod> 1. When Amod = 1, clearly, the modifier has no effect.

Similarly we define the action of the modifier on an inhibitory reaction:

We use a different equation to define steady-state behavior of a system where a single substrate molecule Y is converted into a product Y:

For the rare cases where a non-chemical formulation is needed to describe the system, we provide an alternative algebraic expression for Y∞:

where f is an arbitrary algebraic function and Y1, Y2… are concentrations of other molecules. The use of this algebraic form is discouraged as it weakens the mapping of the model to the underlying chemistry. This form does not admit of modifiers.

Note that these equations are entirely feed-forward: the concentrations of molecules L, M, and Yinput are not affected by their participating in any downstream reactions.

HillTau time course

Equations 1 define the steady state estimate for molecule Y, given a set of molecule concentrations in the preceding layer. We then assume that the approach of the system to steady-state is governed by a simple exponential with characteristic time τ (Fig 1A): where Y(t) is the value of Y at time t. For a simple binding reaction, the time course τ is an experimental observable, and is approximated by τ ~ 1/(kf+kb) where kf and kb are the forward and backward rates for the first-order Hill binding reaction (Figs 1 and S1).

As a slight extension to this, we permit an optional separate time course τ when Y is falling:

The occurrence of different time-courses for buildup and decay is quite common. It happens in a simple binding reaction (S1A, S1C and S1E Fig). It also occurs when there are different chemical steps, such as enzymes with different rates, mediating competing processes for buildup and decay.

HillTau composite form

Putting Eq 1 and Eq 2 together, we have the following closed-form expression for the value of Y at time t + Δt:

The term Y is an optional (positive) baseline level of molecule Y. Δt is the timestep. Note that this is a closed form: Δt can be as large as the end of the simulation.

In the limit of large Δt, we have,

The initial conditions for molecule Y are either specified in the model definition, or as a fallback we estimate the steady-state concentration as in Eq 3.1.

Rate interpretation of HillTau

Formally, HillTau can be seen as an approximation to the following rate equations:

In cases where we have a separate τ2, and Y∞< Y(t):

SBML, and most ODE solvers, will not readily handle this switch between τ and τ2. For the purposes of the equations below we just use τ.

Expanding out Y∞, we have five variants of Eq 4:

Basic activation reaction:

Basic inhibition reaction:

Activation reaction with modifier:

Inhibition reaction with modifier:

Conversion reaction:

Although the HillTau calculations can be done using equations 4.x with a regular ODE solver, the HillTau definition envisions event-driven calculations. Further, the complete HillTau form uses τ2 extensively, which is not handled by the above equations. In principle, the optimization for the HillTau model could just use τ throughout, in which case the reduced HillTau model could be rendered in mass-action form with a reasonable degree of accuracy (S5A Fig).

In Table 1 we summarize the parameters for HillTau. All but the first apply to Reactions.

Simple HillTau models only need species concentrations, and the KA and τ terms for the reactions. The optional terms greatly facilitate the design goals of compactly specifying diverse signaling reactions.

Computing time-evolution and steady-states

To build complex reaction systems, we permit cascading of reactions so that any molecule can be a substrate or equation term in any other reaction. To reiterate, this is a purely feed-forward formulation, so substrates are not affected by any of their downstream targets. We obtain a dependency tree so that on each timestep the updates are carried out in an order which ensures that inputs ripple in order through the cascade of reactions. This may lead to inaccurate estimation of transient responses if the updates are carried out at greater intervals (timestep Δt) than the time-course (τ) of the fastest reaction in the model (S7 Fig). To address this, HillTau assigns an internal timestep Δt which is smaller than the smallest τ in the model. Since one normally performs time-series sampling of reactions at a time finer than the fastest reaction, this restriction usually has little impact on run-time. Further, multi-step systems may include feedback. In such cases the program has to explicitly break the dependency chain. HillTau identifies dependency cycles, picks a reaction based on definition order, and assigns it the next open level in the dependency tree.

A distinct case arises when the HillTau system is used to compute steady state values (e.g., in a dose-response curve). These could ideally be solved by taking an infinitely long time-step. Given the possible presence of feedback, we instead take a long settling time and subdivide it into 10 equal steps so as to allow feedback reactions to also settle.

Motivation for the HillTau formalism

The initial motivation for the HillTau form was the observation that many stimulus-response curves in signaling have a saturating, Hill-like concentration dependence on input strength even if there are multiple intermediate steps [4,44]. Further, many stimulus-response time-courses are visually similar to exponential time-courses. This suggested that a combination of these two properties might be a good approximation even to multi-step signaling cascades.

In order to mathematically support this idea, we considered two of the common motifs in signaling: enzyme-activation of molecules such as phospho-proteins, with a balancing deactivation reaction; and binding of activators to a reagent. Below, we show that the HillTau form achieves a reasonable approximation both to the amplitude and time-course of the response.

First, we considered outcomes of an enzymatic cascade with back reactions. We approximate the rate of production using a Michaelis-Menten- form:

This is balanced by a first-order back reaction:

Then the steady-state at dP/dt = 0 is obtained by combining these:

This is of the same form as Eq 1, showing that a single enzymatic stage in the cascade can be approximated by HillTau. Here we add a subscript 1 to indicate that it is the first reaction in the cascade. Now we stipulate that P1 is the catalyst for the next step, substituting for enzyme E2. This stage results in the formation of product P2:

And so on for multiple steps. Now, suppose that we only have 2 variable inputs to this pathway: the first stage input E and one of the substrates. All other substrates are held fixed. This is a reasonable assumption for HillTau, because any further variable inputs should be treated explicitly either as modulators or as separate reaction steps. Then, all the (Kmn + Sn) terms are constant except the one variable substrate Sv. By combining all the constant terms into Kcascade, we end up with:

This is equivalent to the Hill equation form at the basis of HillTau (See Eq 1.2). We treat inhibition using the same analysis, except resulting in depletion of a substrate (Eq 1.3).

For the time-course, we assume that one of the reactions is rate-limiting. For this step, the rate of formation of product is given by:

This yields an exponential settling curve with a final value P∞., as shown from Eq 5.3. The time-course is given by: where

τ = 1/kr and P0 is the initial value of P.

Note that Eq 5.7 can be rearranged to give Eq 3.

Overall, this approximation yields a consolidated HillTau ‘reaction’ in which we have one stimulus (E, mapping to the activator L in Eq 1), one reactant (Sv mapping to the reagent Yinput in Eq 1), to represent a cascade of enzymatic steps with back reactions.

Next, we consider binding reactions in the pathway. These give the same steady-state Hill Equation form, by definition. From Eq 1.1, setting n = 1, and considering the first reaction generating Y1:

A cascade of similar binding reactions, where Y1 feeds into reaction 2, can also be reduced into the same form: where Rx = R1R2/(R1+KA2) and KAx = KA1.KA2/(R1+KA2)

Here Eq 5.9 has the same form as Eq 1.1 and Eq 5.8.

Using a similar approach, any number of cascading binding reactions will end up fitting the same Hill form. Further, from Eq 5.5 we see that the enzyme/back-reaction steps have a similar Hill-like form. Thus they too can be folded into this cascade.

What is the time-course of this cascaded reaction? As before, we assume that the cascade has one rate-limiting step i. A standard analysis shows that this too has an exponential time-course. which yields the same exponential settling time-course as Eq 5.7: where tau = 1/(L.Kf + Kb) and Y0 is the initial value of Y.

It is important to note that this value of tau has a dependence on a variable, L, and hence the use of a constant value of tau is approximate. This is partially mitigated by utilizing different values of tau for rising and falling phases of the signal Y. During the rising phase, L will have a different (typically larger) mean value than during the falling phase. The use of tau for rising and tau2 for falling phases of the response reflects this.

Thus the steady-state terms for cascading binding and enzymatic reactions can be consolidated into a single HillTau ‘reaction’ step of the Hill form, and the time-course can be approximated by an exponential when there is one rate-limiting step.

Model definition and reference implementation

HillTau reaction systems are set up in a simple JSON format (Fig 1C), for which we have provided a schema. We have implemented a small reference driver program in Python (hillTau.py) that loads the model, runs it with optional stimuli, and plots or saves the output. The hillTau.py file also provides a set of library functions for use in larger programs. An equivalent implementation in C++ using PyBind11 for identical Python bindings is also provided. Three additional utilities are also provided, as described below: for model illustration, model abstraction, and model conversion to SBML. Python scripts for generating the figures in this paper (except Fig 4, which is a schematic only) are also provided in S1 Data and on the GitHub site. Benchmarking was done using the program fig6.py, which calls MOOSE and HillTau through their Python interfaces, and calls COPASI through thePyCoTools Python interface [60]. The output values of multiple benchmarking runs were averaged and used for generating the graphs in fig6_plotting.py. All HillTau and supporting code is licensed under GPL version 3 or later.

Model illustration

We developed a utility htgraph.py, which generates a reaction diagram for HillTau models specified in the.json format. This diagram gives a complete specification of the model topology, in that one can rebuild the structure of the HillTau model by inspection of the reaction diagram, though of course the parameters are not provided in the diagram. htgraph.py uses the dot module of the open-source package graphviz [61] to generate the network graphs.

Model reduction and abstraction: MASH

We provide a utility for performing Model Abstraction from SBML to HillTau: MASH. Briefly, MASH runs the original SBML model with a reference stimulus to explore the key dimensions of its input-output mappings. Typical reference stimuli (built into MASH) include dose-response curves and pulsatile time-series stimuli. MASH then uses standard minimization routines (scipy.optimize library, method “L-BFGS-B”, [62]) to tweak the HillTau model parameters to improve its fit to the original model.

MASH is implemented as a Python script mash.py which uses an ODE model (SBML) as a reference to which it fits a HillTau model. The user provides an initial HillTau model and specifies a series of stimuli to deliver. As part of this the user also defines which are the input molecules, and which are the readouts, and the ODE and HillTau models. MASH generates a topologically identical HillTau model to the original, with parameters optimized to fit. It also reports initial and final scores, expressed as normalized RMS difference between reference model and HillTau. As per Fig 4, the user may introduce additional intermediate steps in the pathway in order to achieve the target model fit. The user specifies a set of stimuli (typically a combination of dose-response and time-series calculations) that explore the model response space. MASH generates a reference response to these stimuli using an ODE solver (MOOSE). The function evaluation for the minimization is carried out by running the HillTau model with the same stimulus, and comparing the HillTau output point-by-point with the reference. The normalized RMS difference over all points is returned as the score. A score of below 0.05 means that on average the original response differs from the HillTau response by less than 5%. MASH uses the scipy.optimize library to tweak the HillTau model parameters to improve the fit, as measured by this RMS score. MASH documentation and examples are provided on the HillTau website. MASH was used to fit the HillTau models for Figs 2, 3 and 5, and the scores are reported. S2, S3 and S4 Figs, and S1 Data specify how these fits were done.

Model conversion to SBML

The utility ht2sbml.py performs a conversion of HillTau models defined in the reference JSON format, into equivalent ODE models defined in SBML. It uses simplesbml (https://github.com/sys-bio/simplesbml for generating the SBML. This uses the forms defined in equations 4.x in Methods. The conversion is approximate on two counts, first, HillTau is an event-driven, not continuous method, and second, HillTau may use different time-courses for rising and falling phases as a reaction proceeds, whereas the ODE form uses only a single time-course. If the HillTau model is generated (for example, after model reduction) such that each reaction only utilizes tau and not tau2, then the approximation is very good. In S5 Fig, we performed ht2sbml.py conversion of 3 HillTau models to SBML and then compared the HillTau output with COPASI calculation of the converted model, under 5 conditions. We obtained an excellent fit (<1% normalized RMS) for an oscillatory model that did not use tau2. The feedback-inhibition model used in Fig 2, which does use tau2, gave a mediocre fit of 29%. The full protein synthesis model was tested under 3 conditions: protein response to BDNF, S6K response to BDNF, and CaMKIII response to Ca. These gave fits of 30%, 26% and 1.7% respectively, though the qualitative response was similar in all cases. Thus the ht2sbml.py conversion works for all HillTau models, but the conversion may be approximate for models which have very different tau and tau2 parameters in their reactions.

HillTau fits to simple mass-action reactions.

Fits are indicated on top of each figure panel. Each of these is a single HillTau ‘reaction’ where ‘input’ is activator in all but Panel E, where ‘input’ is an inhibitor. In all cases the rising phase fits exactly, but in panels A, C and E the falling phase has a different time-course.

(TIF)

Click here for additional data file.

Fit of single-reaction HillTau Model to protein synthesis pathway.

Model is as in Fig 5B. Panels A-F correspond to panels D-I in Fig 5. In all cases protein production rate is readout. Blue plots are reference, orange areHillTau. A: BDNF@3.7 nM + Ca2+@0.2 μM, 900 seconds. B: BDNF@3.7nM, Ca2+@1μM. C: 3 pulses of BDNF@3.7 nM for 5s, coincident with Ca2+@10μM for 1s, pulses separated by 300 s. D: Same as C, but Ca2+ held at baseline of 0.08 μM. E: Dose-response of protein vs. Ca2+, holding BDNF fixed at 3.7 nM. F: Dose-response of protein vs BDNF, holding Ca2+ fixed at 0.08 μM. G: MASH optimization waveform used to fit the HillTau model for protein synthesis.

(TIF)

Click here for additional data file.

Fitting S6K to the protein synthesis pathway model.

HillTau reactions as in Fig 5C. Panels A-F correspond to panels D-I in Fig 5. In all cases activated S6K concentration is readout. Blue plots are reference, orange areHillTau. A: BDNF@3.7 nM + Ca2+@0.2 μM, 900 seconds. B: BDNF@3.7nM, Ca2+@1μM. C: 3 pulses of BDNF@3.7 nM for 5s, coincident with Ca2+@10μM for 1s, pulses separated by 300 s. D: Same as C, but Ca2+ held at baseline of 0.08 μM. E: Dose-response of protein vs. Ca2+, holding BDNF fixed at 3.7 nM. F: Dose-response of protein vs BDNF, holding Ca2+ fixed at 0.08 μM. G: MASH optimization waveform used to fit the HillTau model for S6K activation.

(TIF)

Click here for additional data file.

Fitting CaMKIII to the protein synthesis pathway model.

HillTaureactions as in Fig 5C. Panels A-F correspond to panels D-I in Fig 5. In all cases activated CaMKIII concentration is readout. Blue plots are reference, orange areHillTau. A: BDNF@3.7 nM + Ca2+@0.2 μM, 900 seconds. B: BDNF@3.7nM, Ca2+@1μM. C: 3 pulses of BDNF@3.7 nM for 5s, coincident with Ca2+@10μM for 1s, pulses separated by 300 s. D: Same as C, but Ca2+ held at baseline of 0.08 μM. E: Dose-response of protein vs. Ca2+, holding BDNF fixed at 3.7 nM. F: Dose-response of protein vs BDNF, holding Ca2+ fixed at 0.08 μM. G: MASH optimization waveform used to fit the HillTau model for CaMKIII activation.

(TIF)

Click here for additional data file.

Conversion of HillTau models to SBML, and comparison of the resultant responses simulated in HillTau and COPASI respectively.

A: Oscillator model. This uses only ‘tau’ in its formulation, and fits to within 1%. B. Feedback inhibition model from Fig 2. A 1 uM stimulus is delivered at t = 20, and it lasts till t = 60. This has a mediocre fit of 29%. C-E: Protein synthesis model. C. Comparing protein synthesis response to a BDNF stimulus of 5 nM from t = 2000s to t = 3000s. Fit = 30% is mediocre. D. S6K activation in response to a BDNF stimulus of 5 nM from t = 2000s to t = 3000s. Fit = 26% is mediocre. E. CaMKIII activation in response to a calcium stimulus of 5 uM from t = 2000 to t = 3000s. This fits well, 1.7%. F. HillTau reaction scheme for oscillator model.

(TIF)

Click here for additional data file.

HillTau model schematic for largest model in Fig 6, with 35 HillTau parameters and 11 reactions.

(TIF)

Click here for additional data file.

Dependence of HillTau simulation output on timestep.

In all panels the dashed lines represent the time-series, and the dots represent the sample points for estimating error using the smallest timestep as reference. Accuracy is reported as normalized root-mean square difference from smallest timestep. A: Feedback inhibition. Step stimulus of 1 uM is given at t = 10s, which lasts till t = 50s. 1% accuracy is achieved for dt = 1s. B: feedforward inhibition. Stimulus same as A. 1.5% accuracy at dt = 1s. C: BCM curve. Stimulus of 1 uM is given at t = 10s and stays till the end of the simulation. 1% accuracy at dt = 1s. D: Kholodenko oscillator. Here the system is free-running. 1.2% accuracy at dt = 6s.

(TIF)

Click here for additional data file.

Supplementary code directory.

The zipfile S1 Data with supplementary code expands out into a directory which has Python scripts and model definition files for generating all the simulated components of the figures in the paper, including supplementary figures. Detailed instructions for running the scripts are provided in the README.txt file in the same directory. The scripts should run with Python 3.x, but many figures require additional software installation for the ODE simulators MOOSE and COPASI, as well as some other packages. Details are provided in the README.txt.

(ZIP)

Click here for additional data file.

18 Aug 2021

Dear Prof. Bhalla,

Thank you very much for submitting your manuscript "HillTau: A fast, compact abstraction for model reduction in biochemical signaling networks" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Joanna Jędrzejewska-Szmek, Ph.D.

Guest Editor

PLOS Computational Biology

Kim Blackwell

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Report on the paper "HillTau: A fast, compact abstraction for model reduction in biochemical signaling"

This paper introduces a new methodology and tool for simulating signaling networks. I find the approach an interesting

alternative to mechanistic ODE models. As the author clearly points out, HillTau has less parameters than mechanistic models and

the execution time is several orders of magnitude faster. These crucial advantages are

appropriately exploited by the machine learning features of the present tool.

However, I find the presentation of the tool and methodology rather expeditious. The mathematical foundations are far from

clear and may lead to confusion. The relation with other abstraction methods is also not clearly presented.

Given the added value of the tool and methodology I recommend publication, provided that the following points are clarified:

1) Introduction.

a) The most current approach in machine learning of dynamical signaling networks is based on neural network type models employing

sigmoidal functions (see for instance Nyman et al Plos Comp Biol 2020, PMID: 32667922). A parallel/comparison to these approaches should be drawn.

b) The fast reaction detection invoqued at line 72 is a particular case of more general quasi-equilibrium and quasi-steady state

based model reduction methods reviewed in Radulescu et al Frontiers in Genetics 2012, PMID: 22833754. These reduction methods use graph rewriting and

retain indeed the chemical reaction network formalism as mentioned at line 76.

2) Hilltau schemes.

The figures and examples are useful for understanding what a Hilltau model is, but are not sufficient. I was not able to find in the paper

an algorithmic description about how to generate a Hilltau scheme given a full mechanistic model. Are these abstractions obtained

manually? It would be also interesting to discuss the generation of Hilltau models starting with protein-protein interaction

networks that are available in larger quantity and scale than the mechanistic models.

3) Testing accuracy.

The author has presented the results of extensive benchmarking of his tool. However, I would have liked to see extensive data

not only on the efficiency (Figure 6) but also on the accuracy of HillTau.

4) Domain of validity.

Several properties of mechanistic models are necessarily lost by this approach. For instance, even if multistationarity with

feedback can be reproduced, multistationarity without feedback (resulting from sequestration effects, see for instance Markevich et al J Cell Biol 2004

PMID: 14744999) is lost. There is a bit of confusion here concerning conservation of the steady states. The method conserves steady states

to some extent, but not stoechiometry classes; this subtle difference leads to the lost of multistationarity without feed-back.

It would be useful that the author warns about other properties that may be lost by using

HillTau abstraction.

5) Methods.

a) definitions should be provided for different concepts, for instance input concentration to a signaling step, level of a signaling step, ligand, modifier, etc.

b) It is not clear why Eq 1.7 does not admit modifiers (line 599)

c) rate interpretation of HillTau (lines 637 - 670). I could understand this heuristics, but here one should notice that the ligands and modifiers are considered

constant in the rate interpretation. This assumption may work for small time steps, but not for large time steps. The author seem to care about the size

of the time step only in the presence

of feedback, whereas large time steps may lead to inacurrate dynamics also without feedback. The main issue here is that the dynamics of one variable in a

network is not mono-exponential, but multiexponential. Only when the timescales of the variables are well separated, single variables have mono-exponential

dynamics. This remark can be related to the previous point 4).

d) motivation for the HillTau formalism (lines 688-747). This part of the methods should be entirely revisited as it is neither clear nor rigorous.

Reviewer #2: The author proposes a new coarse grained modeling approach using equations the author calls HillTau which incorporates Hill like equations coupled to an exponential term to model time delays. There could merit in this approach as there is a significant problem in being able coarse grain signaling models and a solution to this problem is needed. I would like to see one thing added to the paper near the beginning which is a step-by-step example of how someone could take, say the MAPK pathway with feedback and convert it into a HillTau model. I note that the author has conveniently provided two python utilities to convert to and from HillTau models and SBML. I wonder if the simulation tooling the author provides is even needed because once in SBML it can be modeling by many other applications, but that is more of an observation than anything else. Finally, the authors states:

“It is also possible to use linear algebraic root-finding to find the steady-state value in one step, but this is prone to numerical challenges even in modern simulators like COPASI”

I don’t believe this. The author cites Clarke 1981 to support this claim but since that publication a lot of work has been done on the problem and so long as the initial starting point is not too far away from the solution, using newton solvers to find the steady state invariably works. NLEQ and Kinsolv are two very robust newton solvers. I agree that problems will arise if there are conserved moieties in the network which are common in signaling pathways. However, this is easily taken care of by doing a row reduction of the stoichiometry matrix first. One can even do a short presimulation (which many of the mainstream simulators do) to give the solvers a good starting point. I am happy to be corrected if the author can find more recent publications that suggest there are problems.

Reviewer #3: This manuscript describes a framework, dubbed HillTau for representing reaction-diffusion systems that is at once both familiar and novel. The framework is familiar in that it is based on Hill equations for steady states and biologically meaningful time constants. The novelty lies in the event-based simulation this approach enables, as the time of the next transition can be calculated using simple arithmetic expressions. This approach vastly reduces the number of parameters relative to mass-action kinetics by avoiding the need to simulate every interaction and instead focus on the relationships between key concentrations. The potential benefit from reducing the number of parameters and avoiding the need for integration through continuous time is vast as this will allow models to incorporate reaction-diffusion based kinetics at scale with minimal computational overhead. The discussion highlights four potential areas of application: model reduction, system abstraction, scaffolds for model fitting, and efficient approximations to complex signaling.

I think this is mostly a very promising paper with an associated well-documented GitHub repository, however I do have two larger concerns: (1) I'm not sure how we're supposed to interpret quality of fit, (2) the order of the paper was confusing with the HillTau model not being explained or motivated until the methods section.

With respect to quality of fit, the main metric used here is RMS error, which gives promising values, but visually the orange (HillTau) and blue (original) look fairly different to me in Figures 3C, 3D, 5D-G, and 5J. In 5G, there is a systematic bias, with the HillTau solution always below the mass-action version. Curiously, there is even visible difference in 2C, which has fairly smooth kinetics. To the extent that the output of these kinetics may feed into some non-linear system, these variations could be significant. The fits may very well be good enough for modeling, but this is not currently self-evident.

As far as the order goes, the results section opens with a brief description of HillTau, but until getting to the methods, there was no concrete example or equation, making it hard to understand what exactly the approach is that's being described and how it allows event-based simulation. The multiple time constants mentioned in line 326 had not been previously mentioned, leading to additional confusion about the nature of the method. All-in-all, the paper does a good job describing the method, but it could benefit from some of that description and motivation happening earlier. Perhaps the HillTau formalization could be thought of as a result, with the mehtods being mostly e.g. how to fit?

Minor:

It is not clear what the numbers in e.g. Figure 2F and 2G (and similarly elsewhere) mean.

Given that a key strength is that the time constants are supposed to be biologically meaningful, it's too bad that the models are being fit using scipy optimization. Is there a principled way to derive the time constants from the mass-action rates?

In at least some of the examples e.g. Examples/PaperFigures/bench_native.py, elapsed time is measured using deltas of time.time(), which gives timings to the nearest 1/64 sec. time.perf_counter() gives higher resolution times that are not susceptible to e.g. changes in system clock. See https://www.webucator.com/article/python-clocks-explained/ for a discussion.

Line 345: "setused" -> "set used" (missing space)

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

16 Sep 2021

Submitted filename: ResponseToReviewers.pdf

Click here for additional data file.

12 Oct 2021

Dear Prof. Bhalla,

Thank you very much for submitting your manuscript "HillTau: A fast, compact abstraction for model reduction in biochemical signaling networks" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Joanna Jędrzejewska-Szmek, Ph.D.

Associate Editor

PLOS Computational Biology

Kim Blackwell

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The revised version of the manuscript addresses all my recommendations and can be accepted.

Reviewer #2: No further comments on the paper.

Reviewer #3: I thank the author for his responses, which largely address my concerns. This manuscript has been strengthened by his edits, especially the addition of a summary of the nature of the HillTau model at the beginning of the results section. In brief, HillTau can be viewed as a phenomenological model that lumps multiple reactions together. Key features are its event-driven nature which allows simulations to make relatively large advances as compared to differential-equation based models, and asymmetric exponential change. Feedback cycles are not directly compatible with the model, but can be approximated by breaking the cycle and using small advances.

Although, honestly, I'm still a little confused because the approach uses an internal time-step so it can't skip from event to event.

The existing implementation on GitHub is pure Python, but the readme promises a faster C++ version that preserves the existing Python API. This will be invaluable.

I still wonder how I would decide if a simplified HillTau model was good enough for my purposes, but the authors points in lines 397-406 (all line numbers here and below are in the version with edits visible) do a good job of addressing the tradeoffs inherent in building "compact models."

Minor:

lines 183-184: A factor of 10 usually gives good convergence. <-- would benefit from knowing what evidence this was based on

MASH is introduced in lines 195-197 , used in 209-210, then seemingly reintroduced in 358-361... would be better if this was consolidated

line 357: readers may not be familiar with FindSim (Viswan et al. 2018); the paper would benefit from a brief description... it's not clear why that approach is best here vs other approaches

lines 457-458: not 100% clear what you're considering state variables here; usually these would be things that have associated rates of change, and then you very much would want to know the initial values. I'm thinking of things like gating variables in Hodgkin-Huxley. There's no initial concentration but there's definitely an initial value that matters.

Line 469: mysterious space in the middle of C++.

Line 914: I think "HillTaul" is a typo, otherwise I missed something major.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: None

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

26 Oct 2021

Submitted filename: Response_v3.pdf

Click here for additional data file.

8 Nov 2021

Dear Prof. Bhalla,

We are pleased to inform you that your manuscript 'HillTau: A fast, compact abstraction for model reduction in biochemical signaling networks' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Joanna Jędrzejewska-Szmek, Ph.D.

Associate Editor

PLOS Computational Biology

Kim Blackwell

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #3: The author has addressed my concerns, and I have no further comments.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #3: No

24 Nov 2021

PCOMPBIOL-D-21-01193R2

HillTau: A fast, compact abstraction for model reduction in biochemical signaling networks

Dear Dr Bhalla,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Katalin Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol