Comparison of rule- and ordinary differential equation-based dynamic model of DARPP-32 signalling network

Previous comparisons

The potential of the RB paradigm has been extensively discussed (Danos, 2007; Chylek et al., 2013, 2014, 2015), exemplified (Danos, 2007) and used to answer novel biological questions across signalling, regulatory and metabolic networks (Wilson-Kanamori et al., 2015; Antunes, Roque & Simoes-de Souza, 2016; Di Camillo et al., 2016; Santibáñez, Garrido & Martin, 2020; Chattaraj, Blinov & Loew, 2021; Nosbisch, Bear & Haugh, 2022, see reviews for earlier models published before 2006—Hlavacek et al., 2006; and from 2007 to 2013—Chylek et al., 2014). These models are often based on chemical reaction networks previously developed for ODEs. Though the RB paradigm was implemented to match the concept of chemical reactions, rules are generators of reaction networks that may lead to diverging results. As the modelling paradigm can affect the underlying model specification, it would be informative to compare simulations of RB and ODE models defined by the same reaction network, originally built to be solved with ODEs. RB and ODE models have been compared on various aspects before. For instance, Blinov et al. (2006) compared network-like model of epidermal growth factor receptor (EGFR) based on reactions first defined in the ODE model of Kholodenko et al. (1999), with simplified pathway-like structure. However, several assumptions underlying the original model were purposefully modified what changed the underlying interaction network, such as the dissociation of EGFR dimers when phosphorylated or bound to other molecules. In addition to contrasting rules to reactions, another intensively studied aspect, finalised with positive conclusions, was whether stochastic simulation can reveal new properties of a system previously modelled as deterministic one (Vlysidis & Kaznessis, 2018; Hahl & Kremling, 2016; Bustos et al., 2018). As simplistic as it may sound, none of the comparisons mentioned were intended to painstakingly disassemble reaction using the rules of a relatively medium-sized model and compare simulation results.

Objectives

This study compares the existing ODE model with the new RB model, where both are based on the same definition of chemical reactions. We specifically ask: (1) how closely can the dynamics of an ODE-based model be replicated with an RB one? (2) If differences between the two are observed, what is the underlying cause? (3) If the dynamics are indeed replicated, what advantages are there to using an RB model? We present the results of the models’ comparison at the level of notation and model dynamics generated under different conditions. The advantages and disadvantages of the two model representations are discussed, alongside suggestions for future research. We hope that this small-scale attempt may nevertheless bring value to the modelling community to better choose between the two different approaches.

DARPP-32 and its interaction network

DARPP-32, Protein Phosphatase 1 Regulatory Inhibitor subunit 1B isoform 1 (PPP1R1B) (NCBI Resource Coordinators, 2016) is an important multistate and intrinsically disordered protein (Marsh et al., 2010) regulating DA-dependent synaptic plasticity in medium spiny projection neurons (MSPN), the predominant striatal cell type, where signalling cascades are activated simultaneously by glutamatergic and dopaminergic stimuli (Beninger & Gerdjikov, 2005). The malfunction of DARPP-32 relates to a number of neurological disorders, e.g., Alzheimer’s disease (Cho et al., 2015), addiction (Philibin et al., 2011), affective disorders (Kunii et al., 2014), and schizophrenia (Kunii et al., 2014; Wang et al., 2015). Glutamatergic and dopaminergic signals are integrated by DARPP-32 which is involved in a complex network of interactions regulated by its multiple phosphorylation sites, of which four are known to have an impact on DARPP-32 itself (Yger & Girault, 2011): Threonine 34 (Thr34), Threonine 75 (Thr75), Serine 137 (Ser137), and Serine 102 (Ser102). The Threonine sites (Thr34, Thr75) have major regulatory roles in signal processing. Whereas the Serine sites (Ser137, Ser102) reinforce Thr34 signal.

Fernandez et al. (2006) studied the integrative effect and sensitivities of DA- and glutamate (Glu)-mediated signals on the DARPP-32 network (Fig. 2A). Their model examined the particular effect of cyclic adenosine monophosphate (cAMP)-pulse followed by Ca ${}^{2+}$ spike trains showing that DARPP-32 is a robust integrator, far more complex than a bistable switch between DA and Glu signals. The authors included two main pathways that mediate these signals, cAMP–PKA–DARPP-32 phosphorylated at Threonine 34 (D34) and Ca ${}^{2+}$–PP2B–DARPP-32 phosphorylated at Threonine 75 (D75). Contrary to the majority of previous models of the Glu and DA signal integration (Lindskog et al., 2006; Gutierrez-Arenas, Eriksson & Kotaleski, 2014; Nair et al., 2016, 2014), DARPP-32 included three phosphorylation sites: Thr34, Thr75 and Ser137. The authors performed two in silico mutagenesis experiments modifying the role of Ser137. The first mutation inhibits site phosphorylation by changing Serine to Alanine at the 137 position (Ser137Ala). The second mutation leads to permanent phosphorylation of the Ser137 site (constSer137).

Model translation

As a single reaction can be written as a single rule, the translation could have been accomplished in a one-to-one manner. However, to fully benefit from rule patterns, multiple reactions can be condensed into fewer rules by removing irrelevant context, i.e., decontextualised. The context of a reaction in the RB model is defined as the information about the agent’s binding sites, partners and internal states. Based on this definition of reaction context, the following criteria guided decisions about condensing reactions into rules. Given a set of reactions of the same type (forward, backward, or catalytic) between the same reactants (agents), if the difference between reactions lies in agent states (internal or binding) that do not change after the transition from reactants to products, and reaction constants (rates) in all these reactions have the same values, then information about agent states does not define reaction conditions; hence, it can be removed from the reaction notation, i.e., a set of reactions become a rule pattern.

Among the least intuitive cases in encoding reactions into rules are complex substrate activations of PKA and PP2B. PKA is activated by the binding of four cAMP molecules, whereas PP2B activation requires four Ca ${}^{2+}$ ions. In other words, multiple molecules of the same type bind substrates on different sites that have to be uniquely named, which requires explicit encoding of all possible binding combinations on four different sides (Danos, Koeppl & Wilson-Kanamori, 2011), called hereafter as combinatorial binding. Translation from deterministic to stochastic rate constants and molecular concentrations to copy-numbers were performed in a standard way (Feret & Krivine, 2012; Sekar & Faeder, 2012). Lastly, the cAMP pulse and the Ca ${}^{2+}$ spiking are reproduced by the addition of molecular copy numbers and modification of rate constants during the simulation, respectively.

Approach to comparison of models

The qualitative comparison of the results of two dynamic models required the simulation of both models in a stochastic scheme and the alignment of the trajectory of the corresponding observables under varying conditions to comprehensively compare the frameworks (Fig. 1).

Rule patterns reduce reaction number in a certain type of model components

Table 1 juxtaposes the total counts of model components in each model. The RB model represents 152 reactions with 132 rules, each parameterised by one of 62 unique rate constants. This number is lower than the total number of rate constants used to parameterise the ODE model (152). The final rule set is more than twice as large as the unique number of rate constants, meaning that more than one rule is parameterised by the same rate constant. The number of molecular species in the RB model, obtained with snapshots capturing the state of the molecular mixture over simulation time (every 10,000th event), is 91 for the competitive RB model, and 137 for the non-competitive one. In both cases, the sum of molecular species is higher than in the ODE model (75). As expected, the number of rules corresponding to reactions is lower, and the number of molecular species is much higher, confirming that expression patterns reduce the number of rules needed to represent a reaction system. Nonetheless, the number of rules is only slightly lower than the number of reactions (152 to 131). If we closely compare models by parts representing more general molecular mechanisms, rule representation reduces the reaction number in some components but extends it in others (Table 2). The reduction occurred only in “DARPP-32 phosphorylation” and “PP2A activation by Ca ${}^{2+}$” components, where combinations of states of DARPP-32 phosphorylation sites do not have to be explicitly named. In contrast, the increase in the number of reactions compared to the reactions occurred in the components “PKA activation” and “PP2B activation”. They both have four sites that bind the same molecules, Ca ${}^{2+}$ and cAMP, respectively.

RB model recapitulates dynamics of ODE model with minor discrepancies

The RB model recapitulates the principal dynamics of the ODE model, albeit there are some observable differences (compare Figs. 3B and 3C). For instance, during the relaxation phase (after the 600th time point), “D34” in the ODE model needs 100 more time steps to reach the second peak, and it is weaker than its RB counterpart (“D34*”). Worth noting is that the standard deviation in the stochastically simulated ODE model reveals a distinctive variation in abundance of the “D34” observable during the relaxation phase. We further use the stochastic trajectories of the ODE model for comparison with the RB model. For a closer examination, traces of 15 observables (defined in Table S1) obtained from ODE and RB simulations were paired and superimposed (Fig. 4). Next to the clear matches (e.g., Figs. 4B, 4E, 4H and 4N), there are discrepancies between paired curves. Five of these 15 observables (Figs. 4C, 4F, 4I, 4J and 4O) are examples of the largest divergence between models following a similar pattern of behaviour. They are directly connected in a chain of activation reactions that begins with Ca ${}^{2+}$ (Fig. 2B). Higher abundance of all Ca ${}^{2+}$ ions present in the system of the ODE model (Fig. 4C) could explain differences between the remaining four observables. However, the trajectory of “all_Ca” remains at the 0 level during steady states rising only in the spiking interval, which resembles the abundance of free Ca ${}^{2+}$ (Fig. 4B). Ca ${}^{2+}$ activates PP2B represented by the trajectory “PP2Bactive”. The higher level of “PP2Bactive” is consistent with the other three observables, suggesting that this is a factor generating the divergences between the models. Based on the curves of the ODE model, we can reason that a stronger activation of PP2B results in proportionally more copies of the unphosphorylated CK1 and phosphorylated D137. This, in turn, increases substrate availability for PP2C; therefore, more copy-numbers of its bound form. This effect is inverted in the trajectories of the RB model.

RB language allows for detailed dissection of observed molecular species

As the “all_Ca” observable trajectory produced by the ODE model is much lower than in the RB model at the steady state, and “PP2Bactive” appears to dictate the higher effect on the other three observables (“D137”, “CK1u”, “_PP2C”), “all_Ca” and “PP2Bactive” are closer analysed in further steps.

According to the reaction system underlying both models, the activated PP2B is a complex of four Ca ${}^{2+}$ ions and PP2B. This detail is not explicitly stated in the variable name of the ODE model. Therefore, to obtain the trajectory of all Ca ${}^{2+}$ ions, the sum of the copy numbers of the molecular species whose variable names contain “Ca ${}^{2+}$”, would have to be replaced by a more thorough analysis of the relevant reaction context of the ODE model. This is not the case in the RB model, where an observable of interest is obtained with an automated procedure that sums the trajectories of molecular species containing the specified expression pattern. Since the trajectory of all Ca ${}^{2+}$ in the RB model includes the ions bound to PP2B, the comparison of “all_Ca” to “all_Ca*” is inaccurate due to a difference in the molecular species included in these observables. A similar inaccuracy, related to the naming of the observables, explains the discrepancy between the time courses of the total number of cAMP observables (Fig. 3). In contrast to the RB model, the multiple copies of cAMP bound to R2C2 are not included in this ODE model time course (Fig. 4A).

To determine the identity of molecular species whose trajectories were summed up in “all_Ca*”, all species containing Ca ${}^{2+}$ in the RB model simulation were isolated from snapshot data ammounting to 24 compared to only 13 in the ODE model (“all_Ca”). These 13 species correspond in molecular composition to 18 of the 24 RB species sampled in total. The six absent species in the “all_Ca” observable are composed of an active form of PP2B containing four Ca ${}^{2+}$ ions, either free or bound to phosphorylated CK1, or DARPP-32 in four different combinations of phosphorylation states. The number of species comprising “all_Ca” observable in the RB model is higher by five because the half-active form of PP2B (bound to two Ca ${}^{2+}$ ions) in the RB model exists in six variants. Whilst in the ODE model, it is represented as a single species, named “PP2BinactiveCa2”. By rerunning the RB model with a new set of observables matching those in the ODE model and superimposing only one of six trajectories of “PP2BinactiveCa2”, we can obtain a close match between the “all_Ca” observables (Fig. S1C). This demonstrates that the differences between the “all_Ca” observables of the two models can be explained by the difference in the number of representations of the molecular species. As the six trajectories have the same dynamics and average levels of abundances, choosing one of them is arbitrary (Fig. S2). Moreover, the distinction between locations of two Ca ${}^{2+}$ ions on the numbered sites of PP2B is irrelevant since all four sites are functionally indistinguishable.

Rate constants of reactions formulated with “combinatorial-binding” notation should be increased to match ODE trajectories

The largest discrepancy between trajectories can be observed in “PP2BinactiveCa2” (Fig. 5M) that for the RB time course was obtained by summing six entities representing a half-active PP2B into one. If divided by six, representing a single variant of half-active PP2B, the trajectory of the RB model becomes lower than the one of the ODE model (Fig. S3M). These six forms of the half-active PP2B suggest that a better fit between the two models can be achieved by decreasing the constant rate of rules that represent the binding of Ca ${}^2+$ to free PP2B. However, the decrease of this rate constant of PP2B intermediate form would lead to a further decrease of copy numbers of other coupled observables that trajectories are lower than in the ODE model in the current parameter setting, e.g., the fully active form of PP2B (Fig. 4J). To further examine this observation, we can return to the comparison of model specifications (Table 2). The four reactions of PP2B activation are represented by 24 rules, explicit in site-specific detail that includes all combinations of Ca ${}^{2+}$ positioning on four sites of PP2B. Moreover, instead of three distinct PP2B species in the reaction representation (inactive, half-active, activated), in the RB model there are eight different molecular species. Thus, the fine-grained representation of species in the molecular mixture slows down the transition from inactive to active PP2B, despite equal presence of the two-step transition encoded in both formalisms. The same can be observed in the second example that required a much larger number of rules, i.e., an activation of PKA. Figure 4D shows that the RB trajectory of the “PKA*” observable also reaches a much lower peak than its ODE counterpart. Accordingly, values of rate constants of rules that the number increased due to the “combinatorial binding” notation in the RB model, i.e., are represented by more variants of species, should be increased to closely match the ones in the ODE model.

RB language facilitates modifications of dynamic models

Site-directed mutations

The ODE model was analysed under two site-directed mutagenesis conditions affecting the Ser137 site. When these perturbations were applied in the RB model, we could observe a close fit in the initial conditions and a general pattern of dynamics based on six key observables (Fig. 6). Thus, RBM allows to emulate experimentally observed perturbations similarly to ODEs.

Effects of the RB framework on the model notation

Encoding reactions into rules slightly reduced the size of model specification and increased counts of molecular species, which confirmed the well-known advantage of rule representation (Hlavacek et al., 2006; Chylek et al., 2014). Closer analysis of reaction subsets representing more general molecular mechanisms showed that reduction in the reaction number is only true for reactions occurring between the same reactants, describe the same transformation, and are parameterised with the same values of rate constants, but differ only concerning binding or internal state of reactants. In this type of reaction, the number of unique reaction rates was equal to the number of rules. An increase in the number of rules representing reactions occurred where the same partner binds an agent at multiple sites. In such reactions, all possible positions and stages of the binding process had to be explicitly encoded. The “combinatorial binding” notation is not a general property of the RB language but characteristic of the Kappa syntax. In the alternative to Kappa RB framework, BNG, a rule is definable with identical sites’ names. This implies that the rule pattern defined for one applies to the others, effectively shortening the rule description (Sekar & Faeder, 2012).

Discrepancies in dynamics between models

As in other similar study (Blinov et al., 2006), the comparison of trajectories showed agreement between model dynamics, with some discrepancies. Their source is not due to the issue of ambiguous molecularity that can arise in Kappa, which requires control on the part of the modeller, as none of the rules in the DARPP-32 model allow the formation of closed rings between agents. The discrepancies mainly appeared due to a lack of precision in the variable names of the ODE model caused by differences in molecular species comprising tracked observables. When the simulation of the RB model was performed with observables exactly matching the ones in the ODE model, almost all paired trajectories fitted perfectly. The only problematic observable involved the “combinatorial binding” notation. Further analysis suggested that activation of proteins encoded with fine detail was slower than in the ODE model. More specifically, although half-active PP2B was more abundant in the RB model, its active form had lower levels compared to the ODE model. This observation may apply to RBM in general because the molecular mixture when simulated with NFsim (Sneddon, Faeder & Emonet, 2011), a variant of the network-free simulator for BNGL, will implicitly contain the same number of molecular species with different binding variants as it is in KaSim. Further, based on the example of antibodies with identical antigen-binding sites, the BNGL developers demonstrated that the equivalent set of rules in both RBM frameworks are parameterised by the same rate constant (Suderman & Hlavacek, 2017). In the BNG framework, however, this rate is automatically scaled by a factor equivalent to the multiplicity of indistinguishable ways to obtain the reaction product (Faeder et al., 2005; Faeder, Blinov & Hlavacek, 2009). Moreover, comparison of the time traces of Kappa- and BNGL-defined rules shows equivalent results under stochastic realisations (Suderman & Hlavacek, 2017). Thus, we argue that “combinatorial binding” expressed as a single rule in BNGL will not change the dynamics observed with Kappa. These specific discrepancies should be taken into account when the reactions and rate constants of RB models are derived from ODE models, and preferably determined experimentally for these reaction types.

A shift in modelling focus with the use of the RB framework

The process of encoding reactions into rules turns attention to questions such as how many binding partners can simultaneously bind a protein. The translation process has shown that information about interfaces of interacting proteins and their alternative states would considerably ease the process of model development by guiding decisions on agents’ signatures. Therefore, data resources that could support RB modelling, such as a source of protein interaction interfaces, post-translational modifications (PTMs), and protein domains. For instance, proteins containing phosphatase catalytic domains are enzymes of dephosphorylation reactions (Sacco et al., 2012). However, such detailed information is not accessible for most molecular agents.

Facilitation of modification and reuse of dynamic models

Molecular species defined in the ODE framework are fixed and definite, whereas, in the RB model, they are a subject of investigation. The determination of created species and their abundances is enabled with snapshots, state vectors of molecular mixture sampled during the simulation. Snapshots have been used as a tool for model exploration in other studies (Sorokina, Sorokin & Armstrong, 2011; Suderman & Deeds, 2013), and their visualisation is included as a standard tool in the Kappa framework (Boutillier et al., 2018). Although, unlike Kappa, BNG allows counting the number of species through network generation, this option is limited by the size of the model. Then, if not sampled, the number of specie types is calculated analytically (Suderman & Deeds, 2013). While this is an approximate approach, by using snapshots, we were able to study the emerging molecular species in detail and precisely track the observables of interest. Furthermore, the automatic identification and merging of time courses into observables makes the RB framework particularly advantageous over the ODE framework, as it avoids errors in the identification of molecules in the modelled system. This is important when the system consists of many molecules and their states are to be analysed in detail. In ODE models, on the other hand, obtaining the details of the molecule composition, hidden in the arbitrary names of the individual molecular species, would require further deconstruction of the reaction system. Though, the precise identification of molecules within species could be performed by parsing the SBML-model encoding, the Fernandez et al. (2006) model web page in the BioModels website (http://www.ebi.ac.uk/biomodels-main/BIOMD0000000153), shows incomplete annotations of molecular species behind variable names, both concerning the actual counts of interactions (e.g., Ca ${}^{2+}$) and their components.

The Fernandez et al. (2006) study does not discuss if DARPP-32 partners bind to the same or different active sites, though DARPP-32 is an intrinsically disordered protein with an unknown binding interface (Dancheck, Nairn & Peti, 2008; Mollica et al., 2016; Engmann et al., 2015; Choy, Page & Peti, 2012). The ODE model specification demonstrates that DARPP-32 forms at most heterodimers. This type of modification in the ODE model would require the arduous extension to much more complex model by enumeration of additional molecular species, the addition of new equations and updating the existing ones. Contrary to this, a definition of such a binding scenario in the RB notation requires the same number of rules, provided that concurrently bound interactors do not influence each other. The RB model was tested with two types of site-directed modifications, demonstrating the framework’s flexibility to reproduce experimentally conducted perturbations. Though the binding site modification effectively changed the model reaction network, it did not affect the model response. Nevertheless, this intervention demonstrated the ease of performing such alterations within the RB framework. Additionally, the pattern notation improves model clarity and provides an intuitive representation of a model akin to a set of chemical reactions rather than equations, potentially improving the learning curve for a modeller-to-be.

Note on simulation time

The ODE simulation remains incomparably faster than the RB simulation. In both deterministic and stochastic settings, COPASI returns ODE results in seconds, while KaSim 3.5 requires ∼40 min¹ to return the RB model results (single core on ThinkPad Core i7-4700MQ; RAM: 16GB). With more recent versions of KaSim (≥4), the execution time could be much shorter thanks to a new solution for the update of observable instances (Boutillier, Ehrhard & Krivine, 2017). Nevertheless, with a sufficiently simple model, network-based simulation can be performed far more efficiently using BNGL, bearing in mind that more complex models will have a higher runtime (Sneddon, Faeder & Emonet, 2011). When generating large-scale models, the type of RB model simulation can be scheduled on a case-by-case basis (Santibáñez, Garrido & Martin, 2020).

Further explorations of the DARPP-32 RB model

There are two main routes for further explore the RB model. The first one is a modification of parameters defining different phases of combinatorially bound Ca ${}^{2+}$ ions to PP2B and cAMP to R2C2. A particular task would be to identify factors by which the binding constants could be modified to counteract the many intermediate variants of these complexes and the lower copies of their activated final forms. Next, it would be interesting to identify conditions under which we could observe a difference in model dynamics after the addition of binding sites to DARPP-32. In the current setup, the lack of difference might be caused by the similarity in occupancy between a single site and all three sites together, as they do not counter each other’s binding properties. The probability of a site being connected depends on copy numbers of reactants and the strength of binding affinities. Reactions in the model are classified as weak, with dissociation rates in the range of $\mu M$. Low-affinity bindings generally lead to lower levels of site occupancy. Moreover, the amount of DARPP-32 molecules exceeds the total counts of all its interactors. Thus, with the current proportions of reactants, all three DARPP-32 sites cannot be saturated. To expose the potential differences in dynamics between the two binding scenarios, the site occupancy could be modified by increasing the size of the reactant pools or reducing the level of DARPP-32. As mentioned by Fernandez et al. (2006), levels of DARPP-32 vary considerably, between $\mu M$ to tens of $\mu Ms$, in the striatum. With the greater availability of single-cell techniques for protein quantification (Lo et al., 2015) it would be worth establishing more precisely the range of DARPP-32 even at the resolution of a dendritic spine (Otmakhov & Lisman, 2012). Estimating variability between cells could also be used to compare the varying levels of phosphorylated DARPP-32 at Thr34 that were observed in the stochastic simulations.

Extending the model to match evolving knowledge on the DARPP-32 interaction network

The main advantage of RBM over ODE is the flexibility to extend and manipulate the model, rather than to rewrite it into a new one. Through a single representation, the community interested in DARPP-32-related molecular systems could pool their knowledge. We believe that the refactored model of the DARPP-32 in a rule-based framework will facilitate the creation of such an updated model in future. It should be noted that the “DARPP-32 events” pathway (R-HSA-180024.3) specified in the REACTOME database is in fact the model of Fernandez et al. (2006) used in this study. So far, only early signalling events of DARPP-32 have been modelled, localised mainly in the cytosol. The Ser102 site regulates nuclear transportation of DARPP-32 representing late signalling events (Stipanovich et al., 2008). The site was omitted in the Fernandez et al. (2006) model due to the lack of evidence that it may be affected by DA or Glu signalling (Girault et al., 1989; Svenningsson et al., 2004). However, a recent study by Nishi et al. (2017) suggests that Glu can decrease the effect of DA signalling by dephosphorylating DARPP-32 at Ser102 causing the accumulation of DARPP-32 in the nucleus, the effect known to be promoted by drugs of abuse (Stipanovich et al., 2008). Moreover, the role of the DARPP-32 interaction network in the nucleus appears to be an important future goal of the recently published ODE model incorporating DARPP-32 (Yapo et al., 2018), indicating the need to explicitly demonstrate the contribution of DARPP-32 to the observed switch-like behaviour of PKA in the nucleus. The inclusion of the Ser102 site seems inevitable to achieve this goal, but extending the current model with two phosphorylation sites may prove to be quite complex. Encouragingly, to calibrate this extended model, measurements experimentally obtained by Nishi et al. (2017) for all four phosphorylation sites could be fed into Pleione (Santibáñez, Garrido & Martin, 2019), a recently developed tool for this task in RBM. Lastly, it would be worth to establish if stochastic simulation can demonstrate new behaviour of this particular molecular system that deterministic simulation has so far not revealed.

Need for formal prioritisation methods of emerging molecular species

RB modelling offers tools for the dissection of emerging molecular species during simulations. As the number of such molecular species increases it becomes difficult to identify which model components are of particular importance to the system. There is great potential to exploit commonly used methods for model exploration, such as sensitivity analysis, to identify critical model features including parameters and model output variables. It would be advantageous to support the modeller’s assumptions with automated methods to prioritise model outputs for downstream analysis and to gain greater insight into the underlying biological systems.