Ceramide-transfer protein-mediated ceramide transfer is a structurally tunable flow-inducing mechanism with structural feed-forward loops.

This paper considers two models of ceramide-transfer protein (CERT)-mediated ceramide transfer at the trans-Golgi network proposed in the literature, short distance shuttle and neck swinging, and seeks structural (parameter-free) features of the two models, which rely exclusively on the peculiar interaction network and not on specific parameter values. In particular, it is shown that both models can be seen as flow-inducing systems, where the flows between pairs of species are tuned by the concentrations of other species, and suitable external inputs can structurally regulate ceramide transfer. In the short distance shuttle model, the amount of transferred ceramide is structurally tuned by active protein kinase D (PKD), both directly and indirectly, in a coherent feed-forward loop motif. In the neck-swinging model, the amount of transferred ceramide is structurally tuned by active PI4KIIIβ, while active PKD has an ambivalent effect, due to the presence of an incoherent feed-forward loop motif that directly inhibits ceramide transfer and indirectly promotes it; the structural role of active PKD is to favour CERT mobility in the cytosol. It is also shown that the influences among key variables often have structurally determined steady-state signs, which can help falsify the models against experimental traces.

Introduction

Located next to the endoplasmatic reticulum (ER) and close to the cell nucleus, the Golgi apparatus [1,2] has a fundamental role in the life of most eukaryotic cells [3-5]: proteins and lipids synthesized at the ER reach the Golgi apparatus, where they are sorted, appropriately modified, labelled for delivery to a specific location and then sent to their target destination. The Golgi apparatus consists of a stack of several cisternae, delimited by lipid membranes, grouped as cis-cisternae (the closest to the nucleus), medial-cisternae (in the central part) and trans-cisternae (the farthest from the nucleus). The outermost cisternae in the cis and trans faces, respectively, constitute the cis-Golgi network (CGN) and the trans-Golgi network (TGN), and handle proteins and lipids that are, respectively, received and released by the apparatus. When a molecule reaches the TGN, it is typically packaged into a vesicle [6,7] that carries it to its destination. However, a molecule can also leave the TGN due to non-vesicular lipid transport mechanisms [8-10], which do not involve vesicle formation. This intracellular lipid transport occurs at the membrane contact sites (MCSs) between the ER membrane and various organelles [11, table 1]: at the MCSs, the external lipid double layers of the two organelles are sufficiently close to allow for the exchange of membrane lipids, which is greatly facilitated by lipid-transfer proteins (LTPs). At the MCSs between the ER and the TGN, the ceramide-transfer protein (CERT) is highly specific for the transport of ceramide (a sphingolipid) [11-15]. In mammalian cells, ceramide is synthesized at the ER and a CERT-mediated transport process brings it to the TGN, where sphingomyelin synthase (SMS) converts it into sphingomyelin: SMS catalyses the conversion of ceramide and phosphatidylcholine (PC) into sphingomyelin (SM) and diacylglycerol (DAG) [16-20].

CERT-mediated non-vesicular transport of ceramide at the ER–TGN MCSs relies on a complex network of regulatory interactions [14,21], whose exact mechanism is still unclear [17,22-24]. Two different models have been proposed: ‘short distance shuttle’ [12,24] and ‘neck swinging’ [3,14,24]. In the short distance shuttle model, the trafficking of ceramide occurs because CERT moves through the cytosol, continuously shuttling between the ER and the TGN membranes [12] (figure 1a). In the neck-swinging model, CERT can extract ceramide from the ER membrane and release it at the TGN membrane only when it is simultaneously bound to both membranes [3,14] (figure 1b). This paper studies the mathematical description first proposed [27] for these two models, which captures the interplay of PKD, PI4KIIIβ and CERT, the key proteins involved in the mechanism, on an average cellular level in mammalian cells. In Weber et al. [27], the two quantitative dynamical models are calibrated by estimating the parameters from experimental data based on Bayesian inference methods, and then validated against experimental data.

Figure 1.

Models for CERT-mediated ceramide transfer at the ER–TGN membrane contact sites: phases of the transfer process in the short distance shuttle model (a) and the neck swinging model (b). Short distance shuttle: CERT is completely unbound and free to move between the two membranes (phase 1); CERT binds to the ER (because its FFAT motif binds to the VAP molecules present in the ER membrane [12,13,22,25,26]) and its START domain binds to ceramide and removes it from the ER (phase 2); CERT is detached from the ER and moves towards the TGN (phase 3); CERT binds to the TGN (because its PH domain binds to the PI4P molecules present on the TGN membrane [13]) and releases ceramide in the TGN (phase 4); the transfer cycle restarts when CERT detaches from the TGN and is again free to move (phase 1). Neck swinging: to transport ceramide, CERT needs to be bound to both the ER and the TGN membranes (phase 1), so that its START domain can bind to ceramide at the ER (phase 2), remove ceramide from the ER and bend towards the TGN (phase 3), and finally release ceramide in the TGN (phase 4).

Here, the two models are analysed and compared with mathematical tools that help better understand the design principles behind CERT-mediated ceramide transfer.

Improving the understanding of this transport phenomenon is of high interest, not only in itself, but also for therapeutic applications. Identifying how to regulate the ER–TGN ceramide transfer (and thus SM synthesis) is fundamental to controlling several related processes and allowing for novel therapies: targeted biotechnological interventions on the TGN regulatory system could, for instance, optimize the secretion rate of therapeutic proteins in mammalian cells [27,28]. As an excessive conversion of ceramide correlates with resistance of cancer cells to drugs and chemotherapy [29], promising anti-oncogenic therapeutic approaches rely on the ability of acting on (and blocking) CERT-mediated ceramide transfer [30]; on the other hand, a lack of CERT can induce cellular dysfunctions [31]. Therefore, to design suitable interventions, the complex interplay of lipids and proteins at the TGN needs to be carefully understood.

Predictive mathematical models are crucial to revealing essential features of biological systems [32]: dynamical systems that describe the functioning of cells and pathways can be analysed with system-theoretic and control-theoretic methods [33-35] to get more insight into the mechanisms of life and assess key properties and behaviours. However, models are inherently plagued by uncertainty, due to the lack of exact quantitative knowledge and the variability of parameters present in nature; nevertheless, biological systems are incredibly robust [36,37] and able to perform their task in the most different conditions [38-40]. Therefore, it is often necessary—and extremely interesting—to assess the behaviour of a biological system based just on limited, qualitative information. A structural (parameter-free) analysis is then performed to find structural properties [32,41-44] that hold regardless of parameter values (within a feasible domain) and exclusively rely on how the key players are interconnected. Structural approaches have been developed to assess whether, for any feasible choice of the parameter values, a biological system enjoys a fundamental property or preserves a fundamental qualitative behaviour (such as, for instance, stability [45-47], oscillations or multistability [48-52], signed input–output influences [53]).

This paper adopts a structural viewpoint to highlight the fundamental differences between the short distance shuttle model and the neck-swinging model of CERT-mediated ceramide transfer, and reveals structural features of the two models that do not depend on parameter values but are rooted in the networked structure of the interactions among key molecules. To this aim, the models are cast into the mathematical framework of flow-inducing networks [54] and a structural algorithm [53] is employed to understand how variations in a variable can modify the steady-state value of other variables, independent of the value of the system parameters.

Models and methods

The models

CERT-mediated non-vesicular ceramide transfer involves the following essential lipids and proteins: the sphingolipid ceramide, which has to be transferred from ER to TGN; the glycerolipid DAG, produced at the TGN [18] and able to indirectly activate the protein PKD [55,56]; the glycerolipid phosphatidylinositol (PI), along with its phosphorylated version phosphatidylinositol 4 phosphate (PI4P), that attracts CERT to the TGN by binding to its PH domain [13] and is thus essential for ceramide transport [57]; protein kinase D (PKD), existing in several isoforms [58,59], two of which, when activated by DAG, are able to phosphorylate PI4KIIIβ and CERT at the TGN [21,55,56,60]; the cytoplasmic lipid kinase phosphatidylinositol-4-kinase-IIIβ (PI4KIIIβ), a protein that, when activated through phosphorylation by PKD, triggers conversion of PI into PI4P, hence attracting CERT to the TGN [61-63]; the cytoplasmic lipid-transfer protein CERT, highly specific for ceramide transfer [11], endowed with multiple functional domains for binding to the ER and the TGN membranes to extract and release ceramide, and with phosphorylation sites that allow the transfer process to be tuned. The PH domain of CERT is associated with binding to the TGN [3,64,65], the FFAT motif with binding to the ER [14,26], the START domain with binding to ceramide [3,66-68], and the SR motif with the regulation of binding affinities [25,60,69].

The considered ordinary differential equation (ODE) models [27] take into account the different possible states of the key molecules PKD, PI4KIIIβ and CERT and the fundamental interactions among them, which are now briefly summarized.

As shown in figure 2a, both PKD and PI4KIIIβ can be present in two forms: PKD, inactive and unphosphorylated (state 1) and PKDpDAG, active and phosphorylated, bound to DAG at the TGN (state 2); PI4KIIIβ, inactive and unphosphorylated (state 3) and PI4KIIIβp, active and phosphorylated (state 4). The lipid PI4P, which recruits CERT to the TGN, is produced by PI4KIIIβp. PKDpDAG, by direct phosphorylation, can either activate PI4KIIIβ or induce the detachment of CERT from the TGN membrane. Also, PKD can be activated by the CERT-mediated ceramide transfer: at the TGN, ceramide and PC are converted into SM and DAG due to SMS [16,18-20], and the increased quantity of DAG promotes the recruitment of PKD at the TGN, hence its activation.

Figure 2.

Key molecules at the interface between the ER and the TGN can be in different states. (a) Inactive PKD (state 1) and active PKDpDAG (state 2); inactive PI4KIIIβ (state 3) and active PI4KIIIβp (state 4). (b) In the short distance shuttle model, CERT can be either bound to ER (state 5), bound to TGN (state 7) or unbound and phosphorylated (state 6). (c) In the neck-swinging model, CERT can be either unbound and unphosphorylated (state 5), unbound and phosphorylated (state 6) or bound to both the ER and the TGN (state 7).

In the short distance shuttle model (figure 2b), CERT can be present in three forms: CERTaER, bound to the ER and unphosphorylated (state 5); CERTp, unbound and phosphorylated (state 6); CERTaTGN, bound to the TGN and unphosphorylated (state 7). CERT is able to deliver a ceramide molecule only when it is bound to the TGN (state 7). After CERT has delivered ceramide to the TGN, it needs to detach from the TGN (via PKD-induced phosphorylation, which converts it to state 6) and then to be dephosphorylated and bound to the ER membrane (state 5) so as to be able to pick up another ceramide molecule at the ER. Then, CERT can detach from the ER and bind again to the TGN (state 7), due to the action of PI4KIIIβp, to deliver the ceramide molecule.

In the neck-swinging model (figure 2c), CERT can be present in three forms: CERTa, unbound and unphosphorylated (state 5); CERTp, unbound and phosphorylated (state 6); CERTaERTGN, bound to both membranes and unphosphorylated (state 7). CERT is only able to transfer ceramide from the ER to the TGN when it is in the double-bound unphosphorylated state 7. PKD-induced phosphorylation detaches CERT from both the TGN and ER membranes (bringing it from state 7 to state 6) and thus interrupts the transfer process. Then, spontaneous dephosphorylation can bring CERT to state 5, while spontaneous phosphorylation brings it back to state 6. The transition from state 6 to state 7 is mediated by the action of PI4KIIIβp, which attracts CERT to the TGN (concurrently, dephosphorylation and binding to the ER also occur), and is required to reactivate ceramide transfer.

The overall process, with the reactions occurring in the two models according to Weber et al. [27], is visualized in the reaction diagrams in figure 3. In the ODE models, the state variables correspond to the concentrations of the species discussed above: x is the concentration of the molecule in state i. The dots indicate time derivatives. Then, the short distance shuttle system is while the neck-swinging system is The positive constants s1, s3 and s5 represent external inflows, while the controlled input signals u1 and u2 are naturally zero and can be artificially set to positive values; all the coefficients represented by Greek letters are positive.

Figure 3.

Reaction diagrams for the short distance shuttle model (2.1) (a) and the neck-swinging model (2.2) (b): x1 is associated with PKD (state 1), x2 with PKDpDAG (state 2), x3 with PI4KIIIβ (state 3), x4 with PI4KIIIβp (state 4); in the short distance shuttle model, x5 with CERTaER (state 5), x6 with CERTp (state 6) and x7 with CERTaTGN (state 7); in the neck-swinging model, x5 with CERTa (state 5), x6 with CERTp (state 6) and x7 with CERTaERTGN (state 7). In both models: an external inflow produces x1, x3 and x5; all species self-degrade; x1 phosphorylates into x2 [at rate x1f1(u1)]; x2 dephosphorylates into x1 [x2f2(u2)]; x4 dephosphorylates into x3 [κx4]; x3 phosphorylates into x4 due to x2 [x3g1(x2)]; x7 phosphorylates into x6 due to x2 [x7g2(x2)]; x6 dephosphorylates into x5 [λx6]. In the short distance shuttle model: x4-produced PI4P converts x5 into x7 [x5g4(x4)]; x1 turns into x2 due to DAG generated by ceramide transfer, which depends on both ceramide-carrying CERT at the ER and on PI4KIIIβp that recruits it at the TGN [x1h(x4, x5)]. In the neck-swinging model: x4-produced PI4P converts x6 into x7 [x5g4(x4)]; x5 phosphorylates into x6 [ϵx5]; x1 turns into x2 due to DAG generated by ceramide transfer, which depends on double-bound CERT [x1g7(x7)].

The quantitative kinetic model built in Weber et al. [27] is based on chemical reaction kinetics, and resorts to linear mass action kinetics for degradation rates, activation of PKD and dephosphorylation of CERT at the ER, and to enzymatic kinetics of the Michaelis–Menten type for regulatory effects across key proteins. In particular, the functions f, g and h are chosen [27] as f(z) = a(1 + bz), g(z) = a(z/(z + b)), h(z, y) = g[yg(z)] = a(ya(z/(z + b))/(b + ya(z/(z + b)))), where a, b, a, b > 0; with this specific choice of the functions, the parameter values are estimated from experimental data, and the resulting models are validated by comparing simulation results with experimental traces.

This paper, conversely, considers a generalization of the models in Weber et al. [27], where now f( · ), g( · ) and h( · , · ) are generic functions that satisfy the following assumption.

Assumption 2.1.

Functions f( · ), g( · ) and h( · , · ) are non-negative and differentiable, with positive partial derivatives in the positive orthant. In particular, f( · ) are strictly positive, also when their argument is zero: f(0) > 0.

This generalization is motivated by the goal of studying and comparing the two models with a structural approach: in the following, a structural mathematical analysis is performed for the two generalized models, which does not rely on specific function or parameter choices and provides parameter-free insight into the biological functioning. Therefore, the results derived in this paper hold for any choice of the positive parameter values and of the model functions f( · ), g( · ) and h( · , · ) satisfying assumption 2.1.

The functional expressions and the parameter values of the original models in Weber et al. [27] are still used in this paper for the numerical simulations reported in §3.5, where the model parameters and the initial conditions used for the simulations are specified. However, no information about the model initial conditions, parameter values and functional expressions is needed for the following theoretical analysis, because the aim is indeed to give results that are completely independent of uncertain model parameters and of the specific choice of the functions (as long as they satisfy assumption 2.1), as well as of the initial conditions: the results always hold if the system structure is preserved. This ensures robustness of the results with respect to all structure-preserving modelling choices.

The framework: flow-inducing networks

Flow-inducing networks are a general class of models introduced in Giordano & Blanchini [54], which are well suited to describe the behaviour of many biochemical systems where a species can boost a reaction without being affected by it (to a first approximation). In mathematical terms, a flow-inducing network is a dynamical system resulting from the interconnection of various dynamical subsystems, or modules. Each module is compartmental, and thus groups positive state variables subject to mass conservation: the sum of the variables in the module is constant, at all time instants. The time evolution of the variables within each module is due to flows among those variables, which can be either spontaneous (proper flows), or tuned by the value of variables belonging to other modules (flows tuned by flow-inducing signals). A minimal example of flow-inducing network is where a (respectively, b, e) is the concentration of the biochemical species A (respectively, B, E). There is a module including A and B, where the sum of the concentrations is constant: . Species E is part of a different module, evolves according to its own dynamics and can influence the A–B module. The individual concentrations of A and B evolve over time due to: the flow from B to A, which is spontaneous, because the rate of conversion of B into A, f(b), depends on the concentration of B only; and the flow from A to B, which is tuned by species E through a flow-inducing signal. In fact, the rate of conversion of A into B, f(a, e), depends not only on the concentration of A but also on the concentration of E, a species that belongs to another module.

A flow-inducing network can be associated with a graph: the graph representing system (2.3) is in figure 4. In the graph representation, the state variables (concentrations of chemical species) are associated with nodes (the blue squares), the flows with arcs (blue pointed arrows) that connect the nodes, and the flow-inducing signals with meta-arcs (red hammer head arrows) that connect a node to an arc. Nodes associated with variables in the same module are grouped by dotted green boxes. Therefore, the blue arcs always connect nodes within a module, inside the dotted green box, while the red meta-arcs always connect a node in a module to an arc in a different module.

Figure 4.

Graph representation of a minimal flow-inducing network, represented by the system (2.3): the state variables are represented by nodes (blue squares), and variables belonging to the same module are grouped by dotted green boxes; the flows among variables within the same module are represented by arcs (blue pointed arrows) that connect the nodes; the flow-inducing signals are represented by meta-arcs (red hammer head arrow) that connect a node in a module to an arc in a different module.

In the models (2.1) and (2.2), there is no mass conservation for the total species. The two considered models of CERT-mediated ceramide transfer exactly fit into the mathematical framework of flow-inducing networks, as shown in Giordano & Blanchini [54], if an external production rate exactly compensates self-degradation for all variables, thus ensuring mass conservation in each of the modules. However, an analysis in the flow-inducing network framework can be performed for the original systems (2.1) and (2.2) even though there is no mass conservation, because this does not alter the core flow-inducing structure of the models (this aspect is discussed in more detail in the electronic supplementary material, §3).

In this paper, the flow-inducing backbone of the models (2.1) and (2.2), captured by the graphs in figure 5, is exploited to get parameter-free insight into the structure of the interactions, so as to understand the design principles underlying the two models and compare their structures.

Figure 5.

Model graphs as flow-inducing networks: (a) short distance shuttle and (b) neck swinging. The main differences between the models are highlighted. In the short distance shuttle model, the CERT-subsystem has a circular structure, because extraction and delivery of ceramide requires the temporary detachment of CERT from the TGN, and both PKDp (state 2) and PI4KIIIβp (state 4) favour this circulation; in the neck-swinging model, it has a cascaded double-switch structure, where the transfer process, related to the amount of CERT in state 7, can be switched off by PKDp (state 2) and switched on by PI4KIIIβp (state 4). Also PKD is activated differently via CERT at the TGN: in the short distance shuttle model, the flow from state 1 to state 2 is tuned by both the amount of ceramide-carrying unphosphorylated CERT at the ER (state 5) and the amount of PI4KIIIβp (state 4), because PI4KIIIβp recruits CERT from the ER to the TGN by producing PI4P; in the neck-swinging model, the flow from state 1 to state 2 is only tuned by CERT in the double-bound state 7. The graphs represent all the interactions in systems (2.1) and (2.2), apart from external inflows s and species self-degradation.

Algorithmic methods: structural steady-state input–output influences

This section surveys structural steady-state influences in dynamical systems; additional details are in the electronic supplementary material, §2. The algorithm proposed in Giordano et al. [53] to efficiently compute structural influences has been used to obtain the structural results in §§3.3 and 3.4.

Consider a generic nonlinear system and with an output y, an input u and an asymptotically stable equilibrium , such that . Then, the steady-state input–output influence [53,70] is the ensuing variation of the steady-state value of the system output, upon a persistent perturbation due to a constant input applied to the system. Its sign can be computed based on the linear approximation of the nonlinear system in a neighbourhood of the equilibrium , with , , : and where J is the Jacobian matrix computed at the equilibrium, while E and H are a column and a row vector that represent, respectively, how the input acts on the system state and how the output depends on the system state in the linearized system. Then, the sign of the steady-state input–output influence [53] turns out to be

A steady-state input–output influence is structurally signed if, upon a perturbation due to a constant input u, the ensuing variation of the steady-state output value has always the same sign as the input (‘+’, structurally positive influence), always the opposite sign (‘−’, structurally negative influence) or is always zero (‘0’, the so-called perfect adaptation [71-73], which is fundamental in many biological systems [74-77]), for any feasible choice of the parameter values; it is undetermined if the behaviour is parameter-dependent (‘?’).

When a persistent additive input is applied to a single state equation and a single state variable is taken as the system output, the results for all the possible input–output pair combinations can be visualized in the structural influence matrix Σ, whose (i, j)th entry Σ expresses the sign of the overall steady-state influence on the ith system variable of an external persistent additive input applied to the dynamic equation of the jth system variable. The entries of matrix Σ can be computed as , where E = E is a column vector whose entries are all zero but the jth, which is 1, and H = H is a row vector whose entries are all zero but the ith, which is 1.

To structurally evaluate the sign of steady-state input–output influences, the vertex algorithm proposed in Giordano et al. [53] can be applied if the system admits the so-called BDC-decomposition [44,45,49,53], namely, its Jacobian can be written in the form , where the matrices M have rank one and the scalars d are positive and related to the system partial derivatives (see [44,53] for details). The outcome of the algorithm is as follows:

— ‘+’ if the influence is structurally positive, for any feasible choice of the parameters;

— ‘0’ if there is structurally perfect adaptation for any feasible choice of the parameters;

— ‘−’ if the influence is structurally negative, for any feasible choice of the parameters;

— ‘?’ if the influence can have a different sign depending on the chosen parameters.

Numerical methods

The computation of the structural steady-state influence with the algorithm discussed in §2.3 has been implemented in Matlab.

Matlab has been used also to numerically simulate the behaviour of systems (2.1) and (2.2), namely, compute the value of the system variables over time, starting from given initial conditions, until a steady state is reached. In particular, the solutions of the systems of ODEs have been computed using the Matlab function ode15s, well suited for stiff differential equations like those in (2.1) and (2.2). To simulate the systems, functions and parameters have been chosen as in Weber et al. [27].

Results

Comments on the models' flow-inducing structure

The two systems (2.1) and (2.2) can be seen as flow-inducing networks, a concept recently introduced in Giordano & Blanchini [54] (see §2.2). The interaction networks for the two models, seen as flow-inducing systems, are visualized in figure 5.

A flow-inducing network is the interconnection of various subsystems, or modules: indeed, both systems (2.1) and (2.2) can be seen as the interplay of three modules, involving variables x1-x2 (PKD-subsystem), x3-x4 (PI4KIIIβ-subsystem) and x5-x6-x7 (CERT-subsystem). Each module is enclosed in a dashed green box in figure 5. In a flow-inducing network, mass flows involve exclusively variables that belong to the same subsystem (cf. the blue arrow head arcs in figure 5); some mass flows among the variables within a subsystem are spontaneous (those from x2 to x1 and from x4 to x3 in both systems (2.1) and (2.2); that from x6 to x5 in the short distance shuttle system (2.1); those from x5 to x6 and from x6 to x5 in the neck-swinging system (2.2)); other mass flows among the variables within a subsystem are tuned by the value of variables that belong to other subsystems, which creates coupling effects among the subsystems (represented by the red hammer head links in figure 5).

In particular, in the short distance shuttle system (2.1): x2, in the first subsystem, induces the flow from x3 to x4 in the second subsystem; x2, in the first subsystem, induces the flow from x7 to x6 in the third subsystem; x4, in the second subsystem, induces the flow from x1 to x2 in the first subsystem; x4, in the second subsystem, induces the flow from x5 to x7 in the third subsystem; and x5, in the third subsystem, induces the flow from x1 to x2 in the first subsystem. Conversely, in the short distance shuttle system (2.2): x2, in the first subsystem, induces the flow from x3 to x4 in the second subsystem; x2, in the first subsystem, induces the flow from x7 to x6 in the third subsystem; x4, in the second subsystem, induces the flow from x6 to x7 in the third subsystem; and x7, in the third subsystem, induces the flow from x1 to x2 in the first subsystem. These flow-inducing effects are a fundamental aspect to capture how CERT-mediated ceramide transfer can be tuned and governed by the key molecules PKD and PI4KIIIβ.

In the PKD-subsystem, the only difference lies in the feedback term from the other subsystems, which in (2.1) depends on x4 and x5, in (2.2) on x7 only, while the PI4KIIIβ-subsystem is described by the same equations in both models. The important difference is in the form of the CERT-subsystem and in the PKD–CERT feedback relation. In both cases, (i) CERT-mediated transport has a positive effect on the activity of PKD; and (ii) active PKD has two main roles: it detaches CERT from the TGN membrane and it activates PI4KIIIβ. Yet, the outcome is completely different.

In (2.1), the CERT-subsystem has a circular structure: extraction and delivery of ceramide requires the temporary detachment of CERT from the TGN. Hence, the larger the circular flow involving CERTaER, CERTp and CERTaTGN, the more ceramide will be transferred. Active PKD favours this ‘circulation’ directly, by steering the transition of CERT from state 7 to state 6, and indirectly, by promoting the activation of PI4KIIIβ, which in turn steers the transition from state 5 to state 7. Therefore, active PKD sustains ceramide transfer in a coherent feed-forward loop motif [32,78,79], and the overall PKD-CERT feedback loop, as suggested in Weber et al. [27], is positive.

In (2.2), the CERT-subsystem has a cascaded double-switch structure: the process can be switched off by active PKD and switched on by active PI4KIIIβ. As only unphosphorylated CERT that is bound to both membranes can transfer ceramide from the ER to the TGN, the larger the concentration of CERT in this state, the more ceramide will be transferred. Now, the overall PKD–CERT feedback loop does not have a neat sign: active PKD switches off ceramide transfer, but it also activates PI4KIIIβ, which switches on ceramide transfer. Hence, the overall net effect depends on the strength and the time constants of the two antagonistic paths, which form an incoherent feed-forward loop [32,78,79].

Positivity, boundedness, equilibria and Jacobian analysis

The systems (2.1) and (2.2) are positive (the variables, species concentrations, cannot become negative) and bounded (the variables cannot diverge), hence they admit an equilibrium point. In fact, for both systems, the trajectories starting from a non-negative initial condition (x(0)≥0, componentwise) lie in the non-negative orthant for all subsequent time instants (x(t)≥0 ∀t); in particular, the trajectories of each subsystem are bounded in a simplex included in the non-negative orthant.

Proposition 3.1.

System (2.1) and system (2.2) are both positive systems. For any non-negative initial condition, their solutions are bounded in the set for fixed positive constants k1, k2 and k3.

Proof.

For both systems (2.1) and (2.2), when x = 0, with , the corresponding derivative is , hence x cannot become negative: x≥0 for all i. Therefore, the non-negative orthant is positively invariant, which proves positivity. To prove boundedness, the sum of the derivatives for each subsystem can be considered (which, interestingly, is the same for both models). For the PKD-subsystem, As s1 is constant, becomes negative when x1 + x2 exceeds a threshold value k1, hence x1 + x2 ≤ k1. Analogously, for the PI4KIIIβ-subsystem, thus x3 + x4 ≤ k2, while for the CERT-subsystem therefore x5 + x6 + x7 ≤ k3. Hence, the state trajectories of each subsystem are bounded in a simplex in the non-negative orthant. ▪

As the solutions are globally asymptotically bounded in a (convex and compact) set , an equilibrium point must exist inside [80-82]. This proves the existence of an equilibrium. Numerical simulations show that, with reasonable choices of the functions and of the parameters (cf. [27]), the equilibrium point is unique and stable, as expected.

The Jacobian matrices associated with the linearization of the two systems around the equilibrium are here reported: the Jacobian JSDS for system (2.1) is in table 1, top, while the Jacobian JNS for system (2.2) is in table 1, bottom. In both cases, in view of assumption 2.1, all of the Greek letters denote positive values.

Table 1.

Jacobian matrix of the short distance shuttle model, JSDS (top), and of the neck-swinging model, JNS (bottom). Monotone subsystems are squared in green (and correspond to the modules inside the dashed green boxes in figure 5), their interactions in red (corresponding to the hammer head red arrows in figure 5). For JSDS, , , , , , , , , , , . For JNS, , , , , , , , , , . All Greek letters denote positive values.

Both the models are the interconnection of monotone subsystems [83-86]: the diagonal blocks in the Jacobian matrix, squared in green in table 1, are Metzler matrices (matrix M is Metzler if it has non-negative off-diagonal entries: [M]≥0 for all i≠j). Each subsystem is also structurally stable, because the corresponding Metzler matrix is column diagonally dominant () and has negative diagonal entries (M < 0,  ∀i). The flow-inducing interactions that constitute the coupling among the subsystems, leading to mass transfer within a subsystem that is tuned by a variable in another subsystem, are squared in red in table 1.

Ceramide-transfer protein-mediated transfer relies on structurally signed input–output influences

Systems (2.1) and (2.2) admit a BDC-decomposition (cf. §2.3 and §2 of the electronic supplementary material): in fact, the Jacobians in table 1 can be written in the form , where d correspond to the Greek letters and M are rank-one matrices. Therefore, the vertex algorithm proposed in Giordano et al. [53] can be applied to structurally evaluate the sign of steady-state input–output influences in these systems and to compute the influence matrices corresponding to the two models, which are

When studying two different models aimed at describing the same phenomenon, it is interesting to compare the outcomes with each other, and also with possible experimental results. Consistency can be investigated as follows: two signed entries are strongly consistent if they have the same sign (namely, both ‘+’, both ‘−’, or both ‘0’); weakly consistent if there is no contradiction because at least one of the two signs is undetermined, ‘?’; inconsistent if there is contradiction (for instance, one sign is ‘+’ and the other is ‘−’, or ‘0’). The following is a formal definition.

Definition 3.2.

Consider two different models of the same phenomenon, model A and model B (or, equivalently, model A and the experimental observations B), both with input i and output o. The two input–output structural influences σ and σ, computed according to the two different models, are strongly consistent if (σ, σ)∈ {( + , + );( − , − );(0, 0)}, weakly consistent if (σ, σ)∈{( + , ?);(?, + );( − , ?);(?, − );(0, ?);(?, 0);(?, ?)}, inconsistent otherwise. The two structural influence matrices Σ and Σ corresponding to model A and model B (assumed to have the same variables representing the same entities) are strongly consistent if all possible pairs of corresponding entries (Σ, Σ) are strongly consistent, weakly consistent if all possible pairs of corresponding entries (Σ, Σ) are either strongly or weakly consistent, and inconsistent otherwise.

The difference between strong and weak consistency is subtle but fundamental. Weak consistency is not inconsistency, but can lead to inconsistency if additional specifications are introduced. While strong consistency is guaranteed for any possible choice of the parameters, hence also for any restriction to a subset of the parameters, weak inconsistency may no longer hold if further information on the models is available and thus the parameter space is restricted. In fact, weak consistency involves either a sign-determined value compared to non-determined value, or two non-determined values. Yet, when restricting to a subset of the parameter space, non-determined values may become determined, and the two models may become inconsistent. Consider, for instance, weak consistency between ‘?’, for one model, and ‘+’, for the other. If additional information on the first model leads to a restriction of the parameter space where the ‘?’ becomes ‘−’, this leads to inconsistency.

The two structural influence matrices ΣSDS and ΣNS are weakly consistent: there is no pair of corresponding entries having inconsistent signs. In particular, 18 out of 49 entries are strongly consistent and represent signed behaviours that are crucial for the system to perform its function: a persistent input that feeds active (respectively, inactive) PKD will increase (respectively, decrease) the concentration of CERT in all its three states, thus revealing that the presence of active PKD is fundamental for ceramide transfer to occur, and will also increase the concentration of PKD, both active and inactive; a persistent input that feeds active (respectively, inactive) PI4KIIIβ will increase (respectively, decrease) the concentration of CERT in all its three states, thus revealing that also the presence of active PI4KIIIβ is fundamental for ceramide transfer, and will also increase (respectively, decrease) the concentration of PI4KIIIβ, both active and inactive.

Also the structural effect on the steady-state value of each variable due to a persistent step input applied through u1 and u2 in systems (2.1) and (2.2) can be computed with the vertex algorithm in Giordano et al. [53]: Expectedly, the signed influence of u2 is always opposite to that of u1. Again, the results are weakly consistent, and the first four of the seven entries are strongly consistent: in fact, in both models it is fundamental that a persistent input u1 always increases the amount of active PKD and active PI4KIIIβ, decreasing the amount of their inactive counterparts.

Ceramide-transfer protein-mediated transfer is a structurally tunable flow-inducing mechanism

The analysis tools provided in Giordano et al. [53] for structural input–output influences can provide more insight into the structural flow-induction mechanism that regulates ceramide transfer. First, it is necessary to define an output that is proportional to the amount of ceramide actually transferred.

For the short distance shuttle model, as the whole cyclic reaction scheme is required to sustain ceramide transfer, the natural output is the flow variable which suitably quantifies the circular flow involving CERTaER, CERTp and CERTaTGN (namely, the sum of the flows involved in the smallest loop that includes TGN-bound CERT). The structural influence on ΦSDS of step-like persistent variations in the input u, i = 1, 2, and of a step-like persistent input affecting the equation of variable x, i = 1, …, 7, computed based on the algorithm in Giordano et al. [53], is summarized as follows (details on the computation are in the electronic supplementary material, §4): The steady-state influence from any input increasing the concentration x2 of active PKD (either due to a persistent transfer u1 or due to a persistent injection directly applied to the equation of x2) to the output ΦSDS is structurally positive: the steady-state flow always increases. Applying a persistent input to the equation of x1 (inactive PKD) yields the same structurally positive influence. Of course, the input u2 has an opposite effect. Then, the system is structurally designed so that the abundance of active PKD regulates the CERT-mediated flow of ceramide. This is consistent with the observation in Weber et al. [27] that PKD is the dominant regulator of CERT-dependent ceramide trafficking.

For the neck-swinging model, conversely, x7 (the concentration of CERT molecules in the double-bound state) is an appropriate output to quantify ceramide transfer, because CERT-mediated ceramide trafficking can occur only when CERT is bound to both the ER and the TGN. A different output could be the flow variable which quantifies the circular flow involving CERTaERTGN and CERTp (namely, analogously to the short distance shuttle case, the sum of the flows involved in the smallest loop that includes TGN-bound CERT). The steady-state input–output influences computed with the algorithm in Giordano et al. [53] can be summarized as follows (details on the computation are in the electronic supplementary material, §4): If x7 is taken as an output, the influence due to a persistent input applied to the equation of variables x3, x4, x5, x6 and x7 is always structurally positive. However, inputs applied to the equation of variables x1 and x2, as well as inputs u1 and u2, yield an undetermined structural influence. So, PI4KIIIβ has a structurally positive effect on ceramide transfer. Conversely, PKD has no structural effect on ceramide transfer in this model.

Conversely, if ΦNS is taken as an output, the influence from PKD is structurally positive: increasing the concentration of active PKD due to a persistent input leads to an increase in the steady-state flow, regardless of the value of the parameters and regardless of the choice of the functions in the model (provided that they satisfy assumption 2.1). Interestingly, in this case the same behaviour occurs if a persistent input is applied not only to the equation of x1, but also to the equation of x2, x3, x4, x5, x6 and x7 (and, therefore, to any group of more variables). The system seems therefore structurally designed to increase mobility of CERT in the cytosol, which is indeed represented by the variable ΦNS.

It seems therefore that also in the neck-swinging model CERT is expected not to remain attached to both membranes in the same location for a long time, but to frequently unbind, travel in the cytosol and then bind again to both the ER and the TGN. Why is this increased mobility needed? The following hypothesis can be suggested: if CERT remains bound to both the ER and the TGN in the same location, the area will be soon depleted of ceramide. Conversely, increasing CERT mobility in the cytosol could help it detach and bind again where there is much more ceramide to be transferred: this would make the transfer mechanism more efficient. This aspect is significant because it suggests that PKD can actually increase ceramide transfer even in the neck-swinging model, because it induces detachment of CERT from the TGN and is thus responsible for its increased mobility.

Remark 3.3.

The outputs in equations (3.3) and (3.4) are a natural choice to represent CERT flow, after embedding of the two models in the flow-inducing framework. Indeed, ΦSDS and ΦNS represent the circular flows that involve TGN-bound CERT (state 7), as shown in the graphs in figures 3 and 5: they are the sum of the flow components associated with the arcs forming the minimal-length loop that involves the variable x7 in the two graphs. The results above point out that, in the short distance shuttle model, ceramide transport is due to circulation (not to the fraction of CERT in the various states, but to its flow between one state and the other), while in the neck-swinging model what mainly counts is the fraction of CERT in state 7, which is the transport-active state.

The structural results point out an interesting difference between the two models. In the short distance shuttle model, PKD is the only structural regulator of CERT-mediated ceramide transfer, and has a structurally positive effect on ceramide flow. Conversely, in the neck-swinging model, a structurally positive regulation is due to PI4KIIIβ; PKD cannot provide a structurally signed regulation, due to the presence of an incoherent feed-forward loop [32] where it directly inhibits ceramide transfer and indirectly promotes it, through PI4KIIIβ. Yet, even in the neck-swinging model, PKD has a structurally positive effect on the mobility of CERT in the cytosol.

Numerical simulations

This section reports simulations of the behaviour of systems (2.1) and (2.2), for functions f, g and h and parameters chosen as in Weber et al. [27] (cf. §2.1). For the sake of completeness, we report here the systems achieved with the choice of functions and parameters as in Weber et al. [27]: the short distance shuttle model becomes where the parameter values are p11 = 0.6467; p12 = 0.0507; p13 = 0.2064; m11 = 2.3412 × 107; s1 = 4.8511 × 105; a11 = 0.9933; a12 = 9.2877; p21 = 8.1888; p22 = 1.8621; m22 = 7191.9; s3 = 2.0894 × 106; a21 = 1.3142; a22 = 0.0030; p31 = 636.5506; p32 = 4.7598; p33 = 12861; m31 = 4.5024 × 109; m33 = 7.6925 × 106; s5 = 139.0702; a31 = 2.2868 × 10−4; a32 = 7.8218 × 10−4; a33 = 3.2251 × 10−4, while the neck-swinging model becomes where the parameter values are p11 = 1.3493; p12 = 7.1440 × 10−4; p13 = 3.1170; m11 = 5.9974 × 105; s1 = 6.2447 × 105; a13 = 1.2879; a14 = 4.2273; a11 = 0.9933; a12 = 9.2877; p21 = 0.5782; p22 = 1.1462; m22 = 2.1550 × 104; s3 = 6.9681 × 105; a21 = 0.0903; a22 = 8.7353; p31 = 3.1330 × 10−4; p32 = 0.0012; p33 = 207.5791; p34 = 516.5629; m31 = 6.5757 × 106; m33 = 7.8586 × 106; s5 = 2122.7; a31 = 0.0046; a32 = 0.0043; a33 = 7.1526 × 10−4.

The evolution of the two systems above, with given initial conditions (reported in the figure captions), is computed until a steady state is reached; then, a persistent positive step input is applied and the system evolves to a new steady state. The sign of the difference between the old and the new steady states is the signed input–output influence.

The sign of the variation in the steady state of x due to an additive step input applied to the equation of variable x is given by the (i, j)th entry of the following matrices: for the short distance shuttle system and for the neck-swinging system. The entries squared in orange are those for which the influence is structurally signed according to the structural influence matrix: for all of these entries, the simulation results are consistent with the expected sign.

The following tables summarize the effect of a persistent step input applied through u1 and u2: for the short distance shuttle system, while for the neck-swinging system Again, all the structurally sign-determined entries, squared in orange, exactly match the expectation. The time evolution of the state variables is shown for the short distance shuttle system in figure 6a, where a persistent step input has been applied through u1 at the time , and for the neck-swinging system in figure 7a, where a persistent step input has been applied through u1 at the time .

Figure 6.

Time evolution of the short distance shuttle system (2.1), when u1 is suddenly increased from 0 to 105 at time . Initial conditions: x1(0) = 5 × 105, x2(0) = 2 × 103, x3(0) = 2 × 106, x4(0) = 9 × 104, x5(0) = 5 × 105, x6(0) = 2 × 103, x7(0) = 2 × 103 (molecules/cell). (a) State variables. (b) Flow variable.

Figure 7.

Time evolution of the neck-swinging system (2.2), when u1 is suddenly increased from 0 to 105 at time t = 750 min. Initial conditions: x1(0) = 5.7 × 105, x2(0) = 6.4 × 103, x3(0) = 2.07 × 106, x4(0) = 5.8 × 104, x5(0) = 4.4 × 105, x6(0) = 1.6 × 103, x7(0) = 6.8 × 104 (molecules/cell). (a) State variables. (b) Flow variable.

Also the influence on the flow variable of step inputs, either applied through u1 and u2 or added to the equation of the variables x, has been computed and is reported for both the short distance shuttle system and the neck-swinging system The input–output influences that are structurally signed according to the theoretical results, squared in orange, are consistent with the expectation. As an example, the time evolution of the flow variable when a step input is applied through u1 is shown in figure 6b for the short distance shuttle system, and in figure 7b for the neck-swinging system.

Concluding discussion

This work has analysed with a structural approach two dynamical systems that describe PKD–PI4KIIIβ–CERT interactions at the trans-Golgi network in mammalian cells, first proposed in Weber et al. [27], which mathematically represent the mechanism of CERT-mediated ceramide transfer as a ‘short distance shuttle’ [12,24] or as a ‘neck-swinging’ process [3,14,24]. The structural analysis has revealed features that inherently distinguish the two models, regardless of the chosen parameter values, and motif-based [32,39] mechanisms that explain the structurally different functioning of the two models.

In particular, in the short distance shuttle model, active PKD is the only structural regulator of CERT-mediated ceramide transfer and its action is enforced through a structural coherent feed-forward loop [32,78,79], which allows PKD to steer ceramide transfer both indirectly, through PI4KIIIβ, and directly. The flow variable that quantifies ceramide transfer is structurally tuned by inputs that affect the concentration of active PKD.

In the neck-swinging model, on the other hand, PI4KIIIβ is the only structural regulator of CERT-mediated ceramide transfer. The effect of PKD on CERT activity is not clear, because PKD is involved in an incoherent feed-forward loop [32,78,79]: it steers ceramide transfer indirectly, by activating PI4KIIIβ, but reduces it directly. The effect of this incoherent feed-forward loop can be seen in the time evolution of the flow variable that quantifies CERT mobility in the cytosol, simulated in figure 7b, which presents a pulse-like trace with a large spike: spiking behaviours are typical in the presence of incoherent feed-forward motifs [32]. Even though CERT-mediated transfer of ceramide is possible exclusively when CERT is bound to both the ER and the TGN, the system seems structurally designed to have a high mobility of CERT in the cytosol, because a persistent input applied to any of the variables induces a structural increase of the associated flow variable. Why? A possible explanation is that, if CERT remains bound to both the ER and the TGN in the same location, the area will be soon depleted of ceramide, while increasing its mobility in the cytosol could help it detach and bind again where there is much more ceramide to be transferred. This suggests that PKD, by inducing detachment of CERT from the TGN, can increase ceramide transfer even in the neck-swinging model.

It is particularly interesting to reveal clear structural features, which are completely independent of parameters, because they are key aspects of a phenomenon that cannot be modified, at least qualitatively: a structurally positive/negative effect will possibly be weaker or stronger depending on the circumstances, but it will never change sign, no matter how the environmental conditions change. If the sign of this induced variation needs to be preserved at all costs, then it must be crucial for the proper system operation: structural effects are pillars in the system architecture. Conversely, parameter-dependent features reveal aspects where flexibility is fundamental: the system must be able to adapt to the circumstances and turn a positive influence into a negative influence, when the environment changes. In this latter case, sensitivity analysis of the model with respect to parameter values becomes crucial to understand how parameter variations affect the qualitative behaviour.

Those proposed in Weber et al. [27] and analysed in this paper are minimal models that capture the difference between two conceptually different mechanisms for CERT-mediated ceramide transfer, short distance shuttle and neck-swinging, by describing the interactions among PKD, PI4KIIIβ and CERT. A parameter-free structural approach that compares the two different mechanisms must rely on simple phenomenological models that capture the essential functioning and allow to draw mathematically grounded conclusions: the models proposed and successfully experimentally validated in Weber et al. [27] are perfect for this type of analysis, and well worth investigating in view of their ability to reproduce experimental data very well, even though many interactions are neglected in the explicit model formulation. Their agreement with data shows that these models are able to implicitly consider many additional effects—not explicitly modelled—through their coefficients, when the parameter values are suitably chosen.

The considered model is focused on the representation of the two short distance shuttle and neck-swinging mechanisms, hence simplified with respect to other aspects. Indeed, many other interactions occur at the ER–TGN MCSs. Another fundamental key player is oxysterol binding protein (OSBP) [87], a protein that, like CERT, is attracted to the TGN by the presence of PI4P and is released from the TGN when phosphorylated by PKD; it is in charge of transporting sterols and also transports PI4P from the ER to the TGN [88-90]. PKD phosphorylation of OSBP has been shown to impair CERT localization at the Golgi [88] and OSBP has been reported to promote CERT translocation to the Golgi complex [60,87]. The apparent coupling between CERT and OSBP lipid transfer reactions at ER–Golgi MCSs [87,91,92] is explained [93,94] on the basis of the OSBP cycle: in fact, as it transports PI4P to the TGN, ‘OSBP could set the tempo for the delivery of precursors of complex lipids to the trans-Golgi, thereby insuring that the concentrations of cholesterol and sphingolipids increase in parallel along the ER–Golgi interface’ [93, p. 841]. ‘The cellular distribution and the ceramide-transfer activity of the protein CERT might depend on OSBP activity’ [94, p. 946]: by controlling the levels of PI4P at the TGN, OSBP seems to control indirectly the localization and the activity of CERT, which needs PI4P to target the TGN. However, ‘further work is now needed to determine the molecular mechanisms governing OSBP auto-inhibition and their relevance for other FFAT-containing LTPs, including CERT’ [94, p. 945]. Also, recent experimental results [95] interestingly show how the flow of sphingolipids, including ceramide, can control the levels and the dynamic turnover of PI4P at the TGN. Therefore, a very important direction for future work is the construction and experimental validation of new models that explicitly include OSBP and PI4P among their variables, to give a broader picture of the interactions ruling ceramide transfer at the trans-Golgi network, beyond the short distance shuttle versus neck-swinging aspect. Given the crucial role of OSBP, this future research direction is particularly relevant: it would be extremely interesting to perform a structural study also of more complex network models also including OSBP.

The ODE models in Weber et al. [27] consider the behaviour on an average cellular level and neglect spatial variations (in fact, the model validation in [27] is based on experiments where total quantities are measured, regardless of the molecule locations in the cell). Partial differential equation models could be built to include the effect of spatial inhomogeneities: modelling choices could lead to sensible differences [96,97].

Although the models analysed here have been validated through extensive comparisons with experimental data, and fit very well the experimental traces [27], the conclusions that can be drawn based on the analysis carried out in this paper could actually be used to try and invalidate (one of) the models from a structural standpoint. In fact, the proposed structural results need to hold always (regardless of how the model parameters are chosen), provided that the model is correct. Therefore, if some experimental results happened to be in contrast with the theoretical predictions that are based on the model structure, this would invalidate the model itself. In particular, structural influences may be used, in combination with experimental data, to falsify a model: if the structural sign of an input–output influence is conflicting with experimental observations, then the model is not suitable to describe the observed phenomenon. This can be of great help for experimental investigators when a model must be assessed. For instance, if the theoretical analysis reveals a structurally positive influence of species A on species B in a model, then—for any choice of the system parameters—if a persistent input is added to species A, the new steady-state value of species B must increase with respect to its previous value. If even a single experiment shows that this does not happen, then the model is invalidated. Notably, the whole model is invalidated, regardless of the parameters. Also, in the case of inconsistency between two models, structural results can help design experiments to discriminate between the two models. Assume that model 1 predicts a structurally positive influence of A on B, while model 2 predicts a structurally negative influence. An experiment showing that the influence can be positive (respectively, negative) will invalidate model 2 (respectively, 1).

This type of structural insight can be useful not only for a better understanding of biological systems and for model falsification, but also for the construction of artificial biomolecular circuits, in synthetic or bottom-up chemical biology. As an example, the design of the biomolecular oscillators [98,99] has been suggested by mathematically grounded structural considerations.

It is worth stressing that, due to the complexity of the models, the large number of involved species and the huge number of interactions among them, it is virtually impossible to see directly the overall influences among the model variables and to infer structural properties simply by inspecting the system, or the associated graph. The proposed algorithmic approach, instead, provides a systematic tool for structural analysis, which enables model falsification or corroborates model validation, and gives additional insight into the system's inherent design.

In this study, all the theoretical predictions on the structural sign of steady-state influences are qualitatively fully consistent both with the numerical simulations of the two models in §3.5 and with the experimental results and extensive simulations reported in Weber et al. [27, fig. 4 and fig. S6]. In fact, fig. 4 and fig. S6 in Weber et al. [27] show the system reactions to perturbations due to u1 or u2, and the qualitative changes in the steady states of the variables shown in the graphs therein are fully consistent with the theoretical structural influences.

The insight into the inherent differences between the two models of CERT-mediated transport at the TGN in mammalian cells, achieved with the proposed structural analysis, will hopefully inspire new experiments that could be decisive in understanding which of the two models better describes the actual mechanism. Also, analysing a suite of models achieved by augmenting or perturbing the existing ones could help combine various features of the models to exactly reproduce the qualitative system behaviour, thus generating a new model of the real process that may be different from all the models that have been already proposed; this is another interesting direction for future work.