Parameter Inference with Bifurcation Diagrams
Year: 2022

Abstract:
Estimation of parameters in differential equation models can be achieved by applying learning algorithms to quantitative time-series data. However, sometimes it is only possible to measure qualitative changes of a system in response to a controlled condition. In dynamical systems theory, such change points are known as bifurcations and lie on a function of the controlled condition called the bifurcation diagram. In this work, we propose a gradient-based approach for inferring the parameters of differential equations that produce a user-specified bifurcation diagram. The cost function contains an error term that is minimal when the model bifurcations match the specified targets and a bifurcation measure which has gradients that push optimisers towards bifurcating parameter regimes. The gradients can be computed without the need to differentiate through the operations of the solver that was used to compute the diagram. We demonstrate parameter inference with minimal models which explore the space of saddle-node and pitchfork diagrams and the genetic toggle switch from synthetic biology. Furthermore, the cost landscape allows us to organise models in terms of topological and geometric equivalence.

heading: 1 Introduction
text: Inverse problems [1] arise in biology and engineering in settings when the model is not fully known and the desire is to match model behaviour to a given set of observations. This helps systematically guide both model and experimental design. While we would like to understand the quantitative details of a system, often only qualitative changes in response to varying experimental conditions can be robustly measured across independent studies [2, 3]. For example, several studies are likely to agree that the human immune system activates above a threshold concentration of a pathogen and deactivates at a lower threshold concentration, but may disagree on the exact quantities of the thresholds or the magnitudes of the immune response. Bifurcation theory provides us a framework for studying these transitions in a manner that is independent of quantitative details [4]. The emerging picture suggests that identification of the qualitative behaviour – the bifurcation diagram – should precede any attempt at inferring other properties of a system [5].
Inferring the parameters of a model directly from a bifurcation diagram is difficult because it is not obvious how multiple parameters in concert control the existence and position of a bifurcation. It could even be impossible for the model to bifurcate in the manner desired. For models with a sufficiently small number of parameters, finding specific bifurcation diagrams is typically done by
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
hand [6]. Several approaches exist to place bifurcations to desired locations once a manifold is present [7–9] yet typically resort to sampling techniques to search for them in the first place [10, 11]. It is always possible to design bespoke goodness of fit measures to find specific model behaviours, for example using the period and phase of limit cycle oscillations [12]. However, this approach does not generalise across a wider set of qualitative behaviours. Progress has been made in cases where model structure and stability conditions are used to refine the search space [13, 14] yet the resulting objectives are still not explicit in the bifurcation targets and also not differentiable. In the emerging field of scientific machine learning [15–17], parameters of structured mechanistic models are favoured over flexible models in larger parameter spaces. A scalable method for navigating the space of bifurcation diagrams would enable design of differential equations with high-level qualitative constraints. Furthermore one could begin organising models according to qualitatively distinct behaviours.
Back-propagation through differential equation solvers has been a breakthrough over the past couple of years [18, 19] that enabled scalable parameter inference for differential equations from trajectory data. Although one could use trajectory data to create the aforementioned qualitative constraints [20, 21] this would entail over-constraining information originating from the kinetics and dynamical transients of the model. Furthermore, such data usually does not contain sufficient information about dynamical transients in order to identify kinetic parameters. Techniques for back-propagating through implicit equation solvers have also been developed [22, 23] although to the best of the authors’ knowledge have not been applied to bifurcation diagrams at the time of writing this paper.
The problem of inferring differential equation parameters against a user-specified bifurcation diagram decomposes into two parts: searching for bifurcating regimes and matching the locations of bifurcation points to desired values. Matching bifurcation locations is a supervised problem where the data are expressed as bifurcations points [8, 11]. Searching for bifurcations is an unsupervised problem because when bifurcations are not present, there is no distance defined between data and prediction [10]. Therefore only properties of the model can be used to start the search. We propose an approach for performing both tasks in an end-to-end fashion. The bifurcation diagram encodes high-level qualitative information defined by state space structures, rather than kinetics. We apply the strategy of implicit layers [22, 23] to calculate gradients. To compute the diagram we use a predictor-corrector method called deflated continuation [24, 25] developed for partial differential equations.
We find that the cost function landscape contains basins that not only allow us to synthesise models with a desired bifurcation diagram but also allow us to organise models in terms of topological and geometric equivalence. We discuss the relevance of this in model selection. In summary, our paper has the following main contributions:
• An end-to-end differentiable method for locating bifurcations in parameter space and then matching their dependency on a control condition to user-specified locations
• Implementation of the method as a Julia package BifurcationInference.jl • Leveraging the cost landscape for a novel way of organising differential equation models in
terms of geometric and topological equivalence

heading: 1.1 Preliminaries
text: Suppose we collected observations along a scalar control condition p ∈ R and conclude that there are specific values of p for which there are qualitative changes in system behaviour. Let D be the set of those values and let us hypothesise that these transitions occur due to bifurcations in the dynamics that drive the underlying mechanism. Let us model the mechanism with a parametrised set of differential equations for states u ∈ RN with a vector function Fθ in a parameter space θ ∈ RM . For the purposes of introducing this work, we will consider the simplest class of bifurcations known as co-dimension one bifurcations not including limit cycles. Therefore D should contain conditions for which we hypothesise changes in multi-stable behaviour. Let the equations be
∂u ∂t = Fθ(u, p) where Fθ : RN+1 → RN (1)
In the context of the differential equations, and not considering limit cycles for now, we show that a static non-degenerate bifurcation can be defined by a set of conditions on the determinant of the Jacobian ∣∣∂Fθ ∂u ∣∣. The determinant of the Jacobian quantifies the rate at which trajectories in a local
patch of state-space u ∈ RN converge or diverge. Let s ∈ R parametrise the curves that trace out the bifurcation diagram. Any location on the curve u(s) and p(s) must satisfy the steady-state of equations (1). Directional derivatives dds along the diagram require the calculation of a vector that is tangent to the diagram (see Supplementary A). The determinant approaching zero along the diagram means that the dynamics of the system are slowing down, which is an important indicator for the onset of a transition between qualitative behaviours. Furthermore, the slowing down must necessarily be followed by a breakdown of stability; for this to be true it is sufficient but not necessary to require that the determinant cross zero with a finite slope, meaning that its directional derivative along the diagram dds ∣∣∂Fθ ∂u
∣∣ is not zero. This is the non-degeneracy condition. The set of predicted values for the control condition P(θ) ⊂ R at which bifurcations occur are defined as
P(θ) := { p | ∃ u : Fθ(u, p) = 0, ∣∣∣∣∂Fθ∂u ∣∣∣∣ = 0, dds ∣∣∣∣∂Fθ∂u ∣∣∣∣ 6= 0} (2)
A proof of how the conditions (2) are necessary and sufficient for static non-degenerate bifurcations is detailed in Supplementary B. The most common bifurcations between steady states, not including limit cycles, are saddle-nodes and pitchforks [26]. Saddle-node bifurcations, which often appear in pairs (Figure 1A) are defined by stable and unstable fixed points meeting and disappearing. Pitchfork bifurcations occur where a single steady state splits into two stable and one unstable steady state (Figure 1B shows an imperfect pitchfork; a perfect pitchfork arises when θ1 = 0). To illustrate these bifurcations, we define minimal models (Figure 1) that span the space of saddle-node and pitchforks, where indeed zero crossings in the determinant with a finite slope define the set of prediction P(θ). The location of these crossings in general may not match the targets D.
For a given set of parameters θ one could compute the set of predicted bifurcations P(θ) using parameter continuation methods [25, 24]. Our goal is to find optimal parameters θ∗ that match predictions P(θ∗) to specified targets D. We must design a suitable cost function L so that θ∗ := argminθL(θ|D) (3) The optimal θ∗ is not expected to always be unique, but is in general a manifold representing the space of qualitatively equivalent models. Ideally, the cost function L should reward θ for which the number of predicted bifurcations is equal to the number of targets, |P(θ)| = |D|. This is especially important in the case where there are no predictions |P(θ)| = 0.

heading: 2 Proposed Method

heading: 2.1 Cost Function
text: To identify parameter sets that give rise to bifurcation diagrams with specified bifurcation points, we propose a cost function that comprises two terms. The role of the error term is simply to reward predicted bifurcations to coincide with the specified target locations. This of course relies on such bifurcations existing. The role of the eigenvalue term is to encourage an optimiser to move towards parameter regimes that do exhibit bifurcations.

heading: 2.1.1 Error term: matching bifurcations to target locations
text: In order for predicted bifurcations p(θ) ∈ P(θ) to match targets p′ ∈ D we need to evaluate an error term |p(θ)−p′|. A naive approach might take an average over the norms for all prediction-target pairs. However this gives rise to unwanted cross-terms and the possibility of multiple predictions matching the same target without any penalty for unmatched targets. Therefore, we choose a geometric mean over the predictions and an arithmetic mean over targets:
E(θ,D) = 1 |D| ∑ p′∈D ∏ p(θ)∈P(θ) |p(θ)− p′| 1 |P| (4)
The error term is only zero when each target is matched by at least one prediction and allows for cases where the number of predictions is greater than or equal to the number of targets |P| ≥ |D|. An alternative approach, which undesirably introduces more hyper-parameters, would be to let each prediction P(θ) represent the centroid of a mixture distribution and use expectation-maximisation to match the centroids to targets D.

heading: 2.1.2 Eigenvalue term: encouraging bifurcations
text: We can see from Figure 1 and definition (2) that predictions p(θ) can be identified by looking for points along the curve where the determinant crosses zero ∣∣∂Fθ ∂u
∣∣ = 0 with a finite slope d ds ∣∣∂Fθ ∂u
∣∣ 6= 0. Using these quantities we can define a positive semi-definite measure ϕθ(s) of zero crossings in the determinant along a curve parametrised by s which we define as
ϕθ(s) :=
( 1 + ∣∣∣∣∣ ∣∣∂Fθ ∂u ∣∣ d ds ∣∣∂Fθ ∂u ∣∣ ∣∣∣∣∣ )−1
(5)
The bifurcation measure ϕθ(s) is maximal at bifurcations and has finite gradients in nonbifurcating regimes (Figure 2). More specifically, the measure ϕθ(s) is one at bifurcation points and goes to zero an odd number of times between bifurcations. This is because ∣∣∂Fθ ∂u
∣∣ must eventually turn around in order to return back to zero, resulting in the directional derivative dds ∣∣∂Fθ ∂u
∣∣ going to zero. Hence the measure ϕθ(s) goes to zero for each turning point (see Figure 2).
On the other hand, as the determinant ∣∣∂Fθ ∂u ∣∣ diverges, we approach regimes far away from any bifurcations and hence ϕθ(s) → 0. Since we would still like to have non-zero gradients with respect to θ in these regimes we designed the measure to go to zero sufficiently slowly.
While the calculation of the determinant is straightforward, its directional derivative requires a tangent vector to the bifurcation curve. Fortunately the tangent vector Tθ(s) at the solution u(s), p(s) anywhere along the curve s can be calculated as the nullspace of the rectangular N × (N + 1) Jacobian
∂Fθ ∂(u, p) ∣∣∣∣ Fθ(u(s),p(s))=0 · Tθ(s) = 0 (6)
This equation guarantees that the tangent vector Tθ(s) is orthogonal to all hyper-planes defined by the components of Fθ. In this setting the dimension of the nullspace is always known, and therefore can reliably be calculated using QR factorisation methods [27].
Equipped with a measure that quantifies the appearance of bifurcations along a bifurcation arc we can define the total measure for a bifurcation diagram as
Ψ(θ) :=
∫ Fθ(u,p)=0
ϕθ(s) ds∫ Fθ(u,p)=0 ds . (7)
Here we denote ∫ Fθ(u,p)=0
ds as the sum of the line integrals in (u, p) ∈ RN+1 defined by the level set Fθ(u, p) = 0 with s being an arbitrary parametrisation of the curves. The total measure Ψ(θ) is normalised such that Ψ(θ) → 1 in the regimes where the controlled condition region p is densely packed with bifurcations. The total measure Ψ(θ) is added to the error term as if it were a likelihood. This defines the cost function as
L(θ|D) := ( |P| − |D| ) log Ψ(θ) + E(θ,D), (8)
The pre-factor |D| − |P| in the eigenvalue term ensures that the gradients are always pushing optimisers towards a state where |D| = |P|. This can be seen as a step-wise annealing of the eigenvalue term until the desired state is reached.

heading: 2.2 Differentiating the cost function
text: To make use of gradient-based optimisers to locate desired bifurcation diagrams, we show here how to differentiate the cost function. First, we note that while individual bifurcations p(θ) depend smoothly on θ, the total number of predictions |P| does not have gradient contributions with respect to θ. Therefore, we can safely drop the dependency in the prediction counter and now proceed in taking gradients with respect to θ knowing that the only dependencies we need to track are for individual bifurcations p(θ) within the definition the error term (4) and the total measure (7). Therefore,
∂L ∂θ = ( |P| − |D| ) λ ∂Ψ ∂θ Ψ(θ)−1 + 1 |D||P| ∑ p′ ∏ p(θ) |p(θ)− p′| 1 |P| ∑ p(θ) ∂p ∂θ (p(θ)− p′)−1 (9)
In a similar vein to back-propagation through neural differential equations [18] we would like to be able to calculate the gradient ∂L∂θ without having to differentiate through the operations of the solver that finds the bifurcation diagram Fθ(u, p) = 0 and the bifurcation locations p(θ). To calculate the gradient of the measure ∂Ψ∂θ we need to differentiate line integrals that depend on θ. Fortunately this can be done by the application of the generalised Leibniz integral rule, details of which can be found in Supplementary C.
The gradient of the bifurcation points ∂p∂θ is found by application of the implicit function theorem to a vector function Gθ : RN+1 → RN+1 whose components represent the two constraints Fθ(u, p) = 0 and ∣∣∂Fθ ∂u ∣∣ = 0. By following a similar strategy to that used by implicit layers [22] we yield an (N + 1)×M Jacobian representing a deformation field [28] for each θ direction. The gradient we are looking for becomes
∂p ∂θ = −p̂ · ∂Gθ ∂(u, p) −1 ∂Gθ ∂θ ∣∣∣∣ Gθ(u,p)=0 where Gθ(u, p) := [ Fθ(u, p)∣∣∂Fθ ∂u ∣∣ ] (10) Here p̂ is a unit vector in (u, p) ∈ RN+1 that picks out the deformations along the p-direction. If we wanted to place the bifurcation at target steady state u′ as well as target control condition p′ we would use the full (N + 1)×M deformation matrix. Calculation of this matrix involves inverting an
(N + 1)× (N + 1) Jacobian ∂Gθ∂(u,p) . Instead of explicitly inverting the Jacobian the corresponding system of linear equations is solved. The determinant of this Jacobian goes to zero in the degenerate case where dds ∣∣∂Fθ ∂u ∣∣ = 0, further justifying our choice of measure Ψ(θ) which discourages the degenerate case.
The cost function is piece-wise smooth and differentiable with undefined gradients only in parameter contours where the number of predictions |P| changes; this is when Ψ(θ) is undefined and the inverse of ∂Gθ∂(u,p) does not exist. Given a set of solutions to Fθ(u, p) = 0 and locations p(θ) the gradient ∂L ∂θ can be evaluated using automatic differentiation methods [29–31] without needing to back-propagate through the solver that obtained the level set Fθ(u, p) = 0 in the forward pass.

heading: 3 Experiments & Results
text: In this section, we apply the method first to minimal examples that can produce saddle-node and pitchfork bifurcations (both N = 1, M = 2), and then a slightly more complex model (N = 2, M = 5) that has multiple parametric regimes producing saddle-node bifurcations. We also demonstrate our method on a model of greater complexity, to convince the reader that the method can be used on more realistic examples with practical significance. In Supplementary D we demonstrate the identification of saddle-node bifurcations and damped oscillations in a model (N = 4, M = 21) of a synthetic gene circuit in E. coli [3].

heading: 3.1 Minimal Models
text: Optimisations of two parameters (θ1, θ2) using simple gradient descent from Flux.jl with learning rate η = 0.01 for the minimal saddle-node and pitchfork models (Figure 1) yield trajectories approach-
ing lines of global minima in the cost function (Figures 3) which represent a set of geometrically equivalent models. Two bifurcation diagrams are geometrically equivalent if the number, type and locations of bifurcations match the specified targets D. We can see that the geometrically equivalent lines are contained within larger basins where the correct number and type of bifurcations are present but do not match the locations of targets D. All models within this basin are in some sense topologically equivalent. This hierarchical classification allows us to identify the set of models that satisfy observed qualitative behaviour [5] before any attempt at inferring kinetic parameters, which is done by choosing a model along the line of geometrically equivalent models.
Optimisation trajectories for the two minimal models appear mostly circumferential. This is because the models were set up such that the radial direction from the origin in θ space mostly scale kinetics whereas the circumferential direction changes the bifurcation topology. This suggests that the gradients of our cost function seek to change model geometry over kinetics.

heading: 3.2 Genetic Toggle Switch
text: In this section we optimise a model where the states share a Hill function relationship with cooperatively n = 2; these models often emerge from mass action kinetics with quasi-steady state approximations and are used to model species concentrations. After re-scaling the equations governing the dynamics of concentrations, the simplified equations for state u1 and u2 become
∂tu1 = a1 + (pu2)
2
1 + (pu2)2 − µ1u1 ∂tu2 =
a2 + (ku1) 2
1 + (ku1)2 − µ2u2 (11)
where ak is the baseline production rate for species k in the absence of the other species. Each species has a finite degradation rate µk. Finally we have two sensitivity constants p and k, one of which is chosen as our control condition. A baseline production rate ak > 1 recovers an inhibitor type hill function for species k and is an activator otherwise. The sensitivities are proportional to the slope of the hill productions. Solving for the steady states, substituting the equation for u1 into u2 and rearranging gives rise to the relationship
k
µ1 =
(1 + ( pµ2u ′)2) √ a2 − u′ (a1 + ( p µ2 u′)2) √ u′ − 1 where u′ := u2µ2 (12)
which reveals that only a1, a2 and the ratio between the sensitivity and degradation parameters, kµ1 , affect the solutions to this equation, and hence the locations of the bifurcations (Figure 4A). In 98% of 800 runs, optimisation using the ADAM optimiser [32] from Flux.jl with learning rate η = 0.1 converged to one of two clusters: mutual activation (a1 < 1, a2 < 1; cluster 1) and mutual inhibition (a1 > 1, a2 > 1; cluster 2) regimes. Example bifurcation diagrams illustrate how the bifurcation curves of each species are positively correlated in mutual activation and negatively correlated for mutual inhibition (Figure 4B).
In order to maintain biological interpretability, optimisation was restricted to the positive parameter regime by transforming the parameters to log-space θ → 10θ. At the beginning of each optimisation run an initial θ was chosen in the log-space by sampling from a multivariate normal distribution with mean zero and standard deviation one.

heading: 3.3 Complexity
text: Performing one iteration of the optimisation requires the computation of the gradient of the cost (9), requiring a computation of the bifurcation diagram with parameter continuation methods, which includes the evaluation of matrix inversions (10). Instead of evaluating the inversions directly, we solve a system of linear equations, applying the same strategy as implicit layers [22, 23]. This leaves us with the computational bottleneck of calculating the determinant of the state space Jacobian, required in both the bifurcation measure (5) and gradient (10). This calculation scales like N2 where N is the number of state space variables (Figure 5A).
For the complexity study, a model was designed so that it is extensible both in the number of parameters M and the number of states N . There are many choices for this; we opted for a model of the form {
∂tu1 = sin 2p− (θ1 sin2p+ 1)u1 ∂tun = un−1 − (µ2n + 1)un 2 ≤ n ≤ N (13)
In this model only the first state u1 defines the shape of the bifurcation diagram, while the remaining states are merely linearly proportional to the first. The parameters µn contain sums of θm allowing us a flexible choice on the number of parameters while maintaining stable solutions for the bifurcation diagram.
While still tractable on laptop computers for states N < 100 our implementation currently does not scale well for partial differential equations where a large the number of states N arises from discretisation of the spatial variables. The only reason we need this determinant is because it
is an indicator of bifurcations. We can address the computational bottleneck by finding a more computationally efficient way of calculating this indicator. One approach would be to take the product of a finite subset of eigenvalues of the system. Note that any more efficient calculation must still permit back-propagation through it.
To demonstrate the benefits of the gradient-based aspect of our method we compare convergence times of gradient descent against a gradient-free approach. We use the Nelder-Mead method from Optim.jl [33] and obtain convergence times as the number of parameters M is increased (Figure 5B). We observe that for our method convergence times scale like √ M compared to M2 for the gradient-free approach.

heading: 4 Conclusion & Broader Impact
text: We proposed a gradient-based approach for inferring the parameters of differential equations that produce a user-specified bifurcation diagram. By applying implicit layers [22, 23] and the generalised Leibniz rule [34] to the geometry of the implicitly defined steady states [35] it is possible to use automatic differentiation methods to efficiently calculate gradients. We defined a bifurcation measure that uses the determinant of the state-space Jacobian as an indicator for bifurcating parameter regimes in the eigenvalue term of the cost function. The gradients of the cost can be efficiently computed using automatic differentiation methods. The computational bottleneck is the evaluation of the state-space Jacobian determinant which limits the implementation to ordinary differential equations.
We demonstrated our approach on models with one bifurcation parameter that can give rise to pitchforks and saddle-nodes. The estimated parameters form distinct clusters, allowing us to organise models in terms of topological and geometric equivalence (Figure 3). In the case of the genetic toggle switch (Figure 4) and a more complex model [3] (Figure D.1) we recovered mutual activation and inhibition regimes. In the more complex model we found a damped oscillatory regime that was not known about in the original paper.
Although we did not consider limit cycles, the bifurcation measure can be extended to detect PoincaréAndronov-Hopf bifurcations alongside changes in stability of fixed points (see Supplementary E for details). This measure enables detection of the onset of damped oscillations and/or the emergence of limit cycles (Figure E.1). Used together with a steady state solver that detects periodic solutions and gradient-based optimisation, we can specify regions of damped oscillation and limit cycles. Our approach generalises naturally to bifurcation manifolds such as limit point curves or surfaces. This is because the normal components of implicit derivatives can still be calculated for under-determined systems of equations [28, 36, 37]. In the case of manifolds it would be more appropriate to use isosurface extraction algorithms rather than continuation to obtain the steady-state manifold. Our approach does not depend on the details of the steady-state solver and therefore can still be applied.
In dynamical systems theory the geometry of state-space determines all of the qualitative behaviours of a system. Our work makes progress towards designing models directly in state-space, rather than the spatial or temporal domain. This is valuable to experimentalists who only have qualitative observations available to them and wish to navigate the space of qualitative behaviours of their system. Our work lies within a trend of progress in the scientific machine learning community, where structured domain-informed models are favoured over flexible models that live in large parameter spaces.

heading: 5 Acknowledgements
text: We would like to acknowledge Kieran Cooney for the fruitful conversations that helped guide the derivations and computational approach. A special thanks go to Romain Veltz and the Julia community for helpful pointers on package development and discussions over Slack. This work was supported by Microsoft Research through its PhD Scholarship Programme and the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, EP/L015854/1).

Review summaries:
In this work, authors suggest a gradient-based semi-supervised approach to infer the parameters of ODEs for producing a bifurcation diagram. A bifurcation measure is defined that uses the determinant of the state-space Jacobian as an indicator for bifurcating parameter regimes in the unsupervised term of the cost function. They demonstrate parameter inference with minimal models which explore the space of saddle-node and pitchfork diagrams and the genetic toggle switch from synthetic biology.

The authors present an approach to inference of parameters of a limited scope of dynamical systems with codimension one bifurcations from user-defined number and location of bifurcation points. The approach solves the parameter inference problem by gradient-based optimization function consisting of two terms. A term matching the provided to the estimated bifurcation locations and a term forcing desired properties of the determinant of the Jacobian of the system. The approach is demonstrated on two simple examples, minimal models and a synthetic toggle switch.