Parameter Inference with Bifurcation Diagrams
Authors: Grisha Szep
Neil Dalchau
Attila Csikász-Nagy
Venue: NeurIPS 2021 Poster (2021)
Keywords: differential equations
bifurcations
parameter inference

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.

Parameter Inference with Bifurcation Diagrams
Gregory Szep King's College London London, WC2R 2LS, UK gregory.o gregory.szepēka. wzep@kcl.ac...
Attila Csikász-Nagy
Pázmány Péter Catholic University Budapest, 1083, Hungary
ceikasznagy@gmaill com
Neil Dalchau
Microsoft Research Cambridge Cambridge, CB1 2FB, UK
ndalchau@gmail. com
Abstract
Estimation of parameters in differential equation models can be achieved by apply- ing learning algorithms to quantitative time-series data. However, sometimes it is only possible to measure qualitative changes of a system in response to controlled condition. In dynamical systems theory, such change points are known as bifurea- tions and lie on function of the controlled condition called the bifurcation diagram. In this work, propose gradient-bar bbasdd approach for inferring the parameters differential equations that produce a user-specified bifurcation diagram. The cost function contains an error term that is minimal when the model bifurcations match specified targets and bifurcation measure which has gradients that push optimisers towards bifurcating parameter regimes. The gradients 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 organise models in terms of topological and geometric equivalence.
be computed
Introduction
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 given set of observations. This helps systematically guide both model and experimental design. While we would like to understand the quantitative details of 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 threshold concentration of 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 framework for studying these transitions in manner that 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 system [5].
Inferring the parameters of model directly from bifurcation diagram is difficult because it is not obvious how multiple parameters in concert control the existence and position a bifurcation. It could even be impossible for the model to bifurcate in the manner desired. For models with sufficienlly small number of parameters, finding specific bifurcation diagrams is typically done by
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

hand [6]. Several approaches exist to place bifurcations to desired locations once 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 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 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- proppappaagttion through differential equation solvers has been 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, 21J 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 mack-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 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 predictor-correctt method called deflated continuation [24, 25] developed for partial differential equations.
We find that the cost function landscape contains basins that not only allow 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 BifuurattionIIfeeencce.
Leveraging the cost landscape for novel way of organising differential equation models in terms of geometric and topological equivalence
1.1 Preliminaries
Suppose we collected observations along a scalar control condition E IR and conclude that there are specific values of p for which there 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 E R with a vector function F' in a parameter space 0 E RRM
For the purposes of introducing this work, we will consider the simplest class of bifurcations known as co- di-mmmsionn one bifurcations not including limit cycles. Therefore 'D should contain conditions for which we hypothesise changes in
behaviour. Let the equations be
Ou Ot
Fo where
Fo N+1 + R
In the context of the differential equations, and not considering limit cycles for now, we show that static non-degenerate bifurcation can be defined by set of conditions on the determinant of the The determinant of the Jacobian quantifies the rate at which trajectories in local
Jacobian OF& du

---

patch of spate-space u E 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 ds along the diagram require the calculation of vector that tangent to the diagram (see Supplementary A). The determinant approaching zero along the diagram means that the dynamies of the system are slowing down, which an important indicator for the onset of a transition between qualitative behaviours. Furthermore, the slowing down must necessarily be followed by a breakdown stability; for this to be true it is sufficient but not necessary to require that the determinant cross zero with finite slope, meaning that its directional derivative along the diagram d aFe is not zero. This is the non-degeneracy condition. The set of predicted values for the control condition P(0) C R at which bifurcations occur are defined as
de '8
aFo P(O) p u: Fo(u, 0, Ou
d OFo
+ 00
ds du
proof of how the conditions (2) are necessary and sufficient for static non-degenerate bifurcations detailed in Supplementary B. The most common bifurcations between steady states, not including limit cycles, are saddle-node 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 IB shows an imperfect pitchfork; a perfect pitchfork arises when 0₁ = 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 finite slope define the of prediction P()). The location of these crossings in general may not match the targets D.
A
B
S
steady states F F
stoondy states
-1-
-2
-2
control condition, p
control condition, p
-1
1 2
-5.0 -2.5 0.0 2.5 5.0
2.5
0.0
OF
du
OF
-25
determinant
determinent
-5
10
-7.5
Figure 1: Illustaaionn of bifurcation diagrams for minimal models of bifurcations. A. Saddle-node bifurcations arise for Fu p) = p + 0 u + 0zu3 when = (g. -1). B. Pitchfork bifurcations arise for (o = 0, + pu + Ozu² when 0 = (1. -1). Targets are illustrated by light yellow vertical lines. Bifurcation curves are shown solid blue and red lines, with lighter shades indicating the determinant crossing zero at locations P()) giving rise to unstable solutions.
For given set of parameters 0 one could compute the set of predicted bifurcations P(0) using parameter continuation methods [25, 24]. Our goal to find optimal parameters 6* that match predictions PP(")) to specified targets D. We must design a suitable cost function L so that
(3) The optimal 0 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(0)] = [D]. This especially important in the case where there are predictions [P()] 0.
0- I
3

---

2 Proposed Method
2.1 Cost Function
To identify parameter sets that give rise to bifurcation diagrams with specified bifurcation points, we propose cost function that comprises two terms. The role of the error term is simply 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.
2.1.1 Error term: matching bifurcations to target locations
In order for predicted bifurcations p(0) E P(0) to match targets p' E we need to evaluate an error term -p'l) A naive approach might take an average over the norms for all prediction-targe 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 geometric mean over the predictions and an arithmetic mean over targets:
E D)
Σ I W'ED P(O)EP(0)
(p()) - p'|
(4)
D
The error term only zero when each target matched by at least one prediction and allows for cases where the number of predictions is greater than equal to the number of targets [P] 2 [D]. alternative approach, which undesirably introduces more byper-parameters, would be let each prediction P()) represent the centroid of mixture distribution and use expectation- cotiatiionmmiiiisttionn to match the centroids to targets D.
2.1.2 Eigenvalue term: encouraging bifurcations
We can see from Figure I and definition (2) that predictions p(0) can identified by looking for points along the curve where the determi- nant crosses zero 8F 8h
0 with finite slope
0F 0. Using these quantities can ds Du
d
define positive semi-definite measure (p((s) of zero crossings in the determinant along a curve parametrised by which define
Measure
(s) == 1+
(5)
Off
The
arelength,
The bifurcation measure (p((s) is maximal at bifurcations and has finite gradients in non- bifurcating regimes (Figure 2). More specifi- cally, the measure (s) is one at bifurcation points and goes zero an odd number of times between bifurcations. This is because s must eventually turn around in order to return back zero, resulting in the directional deriva- tive @P% (s)(s) goes zero for each turning point (see Figure 2).
aF g
0h
determinant
au going zero. Hence the measure
Figure 2: Bifurcation measure (s( and determi- nant iF along the arclength s of two different 2x
On the other hand, as the determinant IFF Puu diverges, we approach regimes far away from any bifurcations and hence (pa(s) + 0. Since we would still like to have non-zero gradients with respect in these regimes we designed the measure to go zero sufficiently slowly.
bifurcation curves demonstrating how maximising the measure along the curve maintains the existing bifurcation marked by circle, while encouraging new bifurcations marked by stars.

---

While the calculation of the determinant is straightforward, its directional derivative requires tangent vector the bifurcation curve. Fortunately the tangent vector (s) at the solution u(s), p(s) anywhere along the curve 8 can be calculated as the nullspace of the rectangular N x (N + 1)
Jacobian
Fo au. p) p(()())
To(s) =
(6)
This equation guarantees that the tangent vector Te(s) is orthogonal to all ppper-planes defined by the components of Fg. In this setting the dimension of the nullspace is always known, and therefore can reliably be calculated using QR factorisation methods [27].
Equipped with measure that quantifies the appearance of bifurcations along a bifurcation are we can define the total measure for bifurcation diagram
(0) Poo(
P ds
Here we denote JF (u ds as the sum of the line integrals in (u. E RN++ defined by the level set p) 0 with being an arbitrary parametrisation of the curves. The total measure (0) is normalised such that (0) + 1 in the regimes where the controlled condition region is densely packed with bifurcations. The total measure (0) is added to the error term as if were likelihood. This defines the cost function as
L = (IP) (Di) log (0 + E(O,
(8)
The pre-factor [D] |P] in the eigenvalue term ensures that the gradients are always pushing optimisers towards state where [D] PPI. This can be seen as a step-wise annealing of the eigenvalue term until the desired state is reached.
2.2 Differentiating the cost function
To make use of gradient- baadeent-based optimisers locate desired bifurcation diagrams, we show here how to differentiate the cost function. First, we note that while individual bifurcations p(0) depend smoothly on 6, the total number of predictions |P] does not have gradient contributions with respect 6. Therefore, we can safely drop the dependency in the prediction counter and now proceed in taking gradients with respect to 0 knowing that the only dependencies we need track are individual bifurcations p(0) within the definition the error term (4) and the total measure (7). Therefore,
aL 00 (IPI (P) X
NV 08
(0)
1
EII
(p()) -p P(O)
Σ
Öp 00 (p(0) -p)- (9)
DIP
In a similar vein to back- propagation through neural differential equations [18] we would like to be at without having to differentiate through the operations of the solver
able to calculate the gradient 20 that finds the bifurcation diagram Fø(u, p) = 0 and the bifurcation locations p(0). To calculate the 3ø we need to differentiate line integrals that depend on A. Fortunately this
gradient of the measure a+
can be done by the application the generalised Leibniz integral rule, details of which can be found in Supplementary C.
Op is found by application of the implicit function theorem to The gradient of the bifurcation points 50 vector function G. N+1 - N+1 whose components represent the two constraints (u. p) = 0 and
& The 0. By following a similar strategy to that used by implicit layers [22] we yield an
(N + 1) x .M Jacobian representing a deformation field [28] for each 0 direction. The gradient we
are looking for becomes
dp
aG.
aGo
where Go(u,)) = Fo p)]
- Hu,
(10)
00
00
Go(u,p)=0
au
Here is a unit vector in (u,p) E N+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) x M deformation matrix. Calculation of this matrix involves inverting an
5

---

(N + 1) x (N + 1) Jacobian aGe Instead of explicitly inverting the Jacobian the corresponding 0(u.p)"
system of linear equations is solved. The determinant of this Jacobian goes to zero in the degenerate case where d 0/ as Du
O, further justifying our choice of measure (0( which discourages the
degenerate case.
The cost function is piece-wiss smooth and differentiable with undefined gradients only in parameter contours where the number of predictions [P] changes; this is when (0) undefined and the inverse 8Gg does not exist. Given
solutions Fo (u, p) and locations p(0) the gradient aL can be evaluated using automatic differentiation methods [29-31] without needing to Back-propagatt through the solver that obtained the level set (g(u. 0 in the forward pass.
of 8(..p)
3 Experiments & Results
In this section, we apply the method first minimal examples that can produce saddle-node and pitchfork bifurcations (both N = 1, M = 2), and then slightly more complex model (N = 2, M = 5) that has multiple parametric regimes producing saddle-node bifurcations. We also demonstrate our method on model of greater complexity, to convince the reader that the method can be used on more realistic examples with practical significance. Supplementary D we demonstrate the identification of saddle-node node bifurcations and damped oscillations in model (N 4, M 21) of a synthetic gene circuit E. coli [3].
0/
slondly stendyydtte states Pi(u, 20
Cost L(OD 92)
parometer,
82
steadystales Palts M
Coss L (1)
-10
5.0
pammnnter,
Figure 3: Saddle-node Fe(u. p) p and pitchfork Fu p) 01 - up 2zu3 optimised with respect to so that predicted bifurcations P(0) match targets D in control condition The right panel shows bifurcations diagrams for the three optimal 0* marked by stars the left panel. The optimisation trajectories in white follow the gradient of the cost, approaching the black lines of global minima in the left panel
3.1 Minimal Models
Optimisations two parameters (0,, 02) using simple gradient descent from Flux. ju with learning rate 77 = 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 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 that the geometrically equivalent lines 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 inferring kinetic parameters, which 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 0 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.
3.2 Genetic Toggle Switch
In this section optimise model where the states share a Hill function relationship with co- operatively 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 govern- ing the dynamics of concentrations, the simplified equations for state and M2 become
+ (puz)² 0u + (pug)² 14111
+ (ku)) 422 1 (kuz)²
0uu2
(11)
where ag is the baseline production rate for species is in the absence of the other species. Fach species has a finite degradation rate f'g- Finally we have two sensitivity constants and k; one of which is chosen
control condition. A baseline production rate 0g 1 recovers 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 ij into u2 and rearranging gives rise to the relationship
k
(1 (1+ (2 w
where u' := 42
(12)
M1
(ar + (
Parameter Estimates 0
B
10
1
steady statos
Sensitivitty Degredation Ratio kip,
10
2
4
6
10
2
16
stoudy states
10"2
10
2 control condition,
4 6
8 10
Basal Production Myy G
Figure 4: Bifurcation inference for the two-state model (11). A. Optimal parameter estimates 6* for the targets DD {4, 5) reveal two clusters of qualitatively different regimes: mutual activation
(a) < 1; cluster 1) and mutual inhibition (a) > 1; cluster 2). B. Example bifurcation diagrams indicate positively and negatively correlated dependencies between the two model states, as function of the control condition.

---

which reveals that only 01, and the ratio between the sensitivity and degradation parameters, k 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. j1 with learning rate 1j = converged one of two clusters: mutual activation (a1 1, d2 < 1; cluster 1) and mutual inhibition (a) > 1, 02 > 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 0 + 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.
3.3 Complexity
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 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 N where N is the number of state space variables (Figure 5A).
A
B 10
y! Cost Gradient
Gradient
Neider-Mead
Convergence Time sec.
Exenution Time sec
10
20
40
States N
P Parameters M
Figure A. Execution time (time to calculate cost gradient) with respect to states N. Convergence times (the time it takes to find and match a bifurcation to within 1% of a specified target) with respect to the number of parameters M, comparing against gradient-fier approach: Nelder-Mead Calculations were performed on an Intel Core 17-6700HQ CPU @ 2.60GHz x 8 without GPU acceleration.
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 many choices for this; we opted for model of the form
{
0.u sin 2p (01 sin + 1)u1 a,u, 1n (4,) 1)un
(13)
2 25n5N
In this model only the first state 2y defines the shape of the bifurcation diagram, while the remaining states are merely linearly proportional to the first. The parameters μ4, contain sums of 6., allowing us flexible choice 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 large the number of states N arises from discretisation of the spatial variables. The only reason 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 take the product of a finite subset of eigenvalues of the system. Note that any more efficient calculation must still permit haok-propagatton through it.
To demonstrate the benefits of the gradient-based 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 j1 [33] and obtain convergence times as the number of parameters M is increased (Figure SB). We observe that for our method convergence times scale like v M compared to M² for the gradient-frre approach.
Conclusion & Broader Impact
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 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.I) 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 F. 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 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 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 trend of progress in the scientific machine learning community, where structured comain-informed models are favoured over flexible models that live in large parameter spaces.
5
Acknowledgements
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/LO15854//

---

References
[1] U. G. Abdulla, R. Poteau, A. Binder, H. W. Engl. C. Flamm, P. K. A. Ugler. J. Lu, S. M. A. Uller, and P. Schuster, "Inverse problems in systems biology," Inverse Problems, vol. 25, p. 51, 2009.
[2] J. Tyson, K. Chen, and B. Novak, "Network dynamics and cell physiology," Nature reviews Molecular cell biology, vol. 2, no. 12, pp. 908-916, 2001.
[3] Grant, G. Szep, O. Patange. J. Halatck, V. Coppard, A. Csikász-Nagy. J. Hascloff, J. Locke, N. Dalchau, and A. Phillipss "Interpretation of morphogen gradients by synthetie bistable circuit," Nature Communications, vol. 11, no. 1, 2020.
[4] Y. A. Kuznetsov, "Topological Equivalence, Bifurcations, and Structural Stability of Dynamical Systems," in Elements of Applied Bifurcation Theory. pp. 39-76, Springer New York, 2004.
[5] M. P. H. Stumpf and E. Roesch, "Parameter inference dynamical systems with co-dimension bifurcations," Royal Society, vol. 6, no. 10, 2019.
[6] A. Csikász-Nagy, D. Battogtokh, K. C. Chen, B. Novák, and J. J. Tyson, "Analysis of a Generic Model of Eukaryotic Cell-Cycle Regulation," Biophysical Journal, vol. 90, pp. 4361-4379, 2006.
[7] K. Iwasaki and Y. Kamimura, "An inverse bifurcation problem and an integral equation of the Abel type," Inverse Problems, vol. 13, pp. 1015-1031, 1997.
[8] J. Lu, H. W. Engl, and P. Schuster, "Inverse bifurcation analysis: Application simple gene systems," Algorithms for Molecular Biology, vol. 1, pp. 1-16, 2006.
[9] I. Dobson, "Distance to Bifurcation in Multidimensional Parameter Space: Margin Sensitivity and Closest Bifureations," in Bifurcation Control: Theory and Applications, pp. 49-66, Springer, Berlin, Heidelberg, 4 2004.
[10] V. Chickarmane, S. R. Paladugu, F. Bergmann, and H. M. Sauro, "Bifurcation discovery tool," BIOINFORMATICS APPLICATIONS NOTE, vol. 21, no. 18, pp. 3688-3690, 2005.
[11] E. D. Conrad, J. Tyson, R. Laubenbacher, J. Phillips, and M. Renardy, "Bifurcation Analysis and Qualitative Optimization of Models in Molecular Cell Biology with Applications the Circadian Clock," Virginia Tech, 4 2006.
[12] J. C. Locke, A. J. Millar, and M. S. Turner, "Modelling genetic networks with noisy and varied experimental data: the circadian clock in Arabidopsis thaliana." Journal of theoretical biology, vol. 234, no. 3, pp. 383-393, 2005.
[13] Otero-Muras and J. R. Banga, "Optimization based prediction of fold bifurcations nonlinear ODE models," IFAC-PupersOnLiie,, vol. 51, pp. 485 490, 2018.
[14] I. Otero-Muras, P. Yordanov, and J. Stelling, "A method for inverse bifurcation of biochemical switches: inferring parameters from dose response curves," BMC Systems Biology, vol. 8, p. 1 114, 2014.
[15] C. Rackauckas and Q Nie, "Differentialequaaions.-- performant and feature-rich ecosystem for solving differential equations in julia," Journal of Open Research Software, vol. 5, no. 1, 2017.
[16] C. Rackauckas, Y. Ma. V. Dixit, X. Guo, M. Innes, J. Revels, J. Nyberg, and V. Ivaturi, "A comparison of automatic differentiation and continuous sensitivity analysis for derivatives of differential equation solutions," arXiv preprint arXiv:1812.01892. 2018.
[17] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, and A. Ramadhan, "Universal Differential Equations for Scientific Machine Learning." arXiv, 2020.
[18] R. Q. Chen, Rubanova, J. Bettencourt, and D. Duvenaud, "Neural Ordinary Differential Equations," NIPs, vol. 109, pp. 31-60, 2018.
10

---

[19] C. Rackauckas, M. Innes, Y. Ma, Bettencourt, L. White, and V. Dixit, "DiffEgFFuuxii - A Julia Library for Neural Differential Equations," ArXiv, 2019.
[20] Ranciatii C. Viroli, and Wit, "Bayesian Smooth-and-Match strategy for ordinary differential equations models that are linear in the parameters." ArXiv, 2017.
[21] F. Khadivar, I. Lauzana, and A. Billard, "Learning dynamical systems with bifurcations," Robotics and Autonomous Systems, vol. 136, p. 103700, 2 2021.
[22] A. Look, S. Doneva, M. Kandemir, R. Gemulla, and J. Peters, "Differentiable Implicit Layers," ArvXiv, 10 2020.
[23] S. Bai, J. Z. Kolter, and V. Koltun, "Deep Equilibrium Models," arXiv, 9 2019.
[24] E. Farrell, C. H. L. Beentjes, and A. Birkisson, "The computation of disconnected bifurcation diagrams," arXiv, 2016.
[25] R. Veltz, "BifuucattinnKitjl"" tech. rep., Inria Sophia-Antipolis 7 2020.
[26] M. Haragus and G. looss, Local Bifurcations, Center Manifolds, and Normal Forms in Infinite- Dimensional Dynamical Systems. London: Springer London, 2011.
[27] Z. Drmač and Z. Bujanović, "On the Failure of Rank-Revealing QR Factorization Software A Case Study." ACM Transactions on Mathematical Software, vol. 35, 2008.
[28] Jos and Schmidt, "On the velocity of an implicit surface," ACM Transactions on Graphics vol. 30, pp.
[29] J. Revels, M. Lubin, and T. Papamarkou, "Forward-Mode Automatic Differentiation in Julia," ArvXiv, 2016.
[30] M. Innes, E. Saba, K. Fischer, D. Gandhi, M. C. Rudilosso, N. M. Joy, T. Karmali, A. Pal, and V. Shah, "Fashionable Modelling with Flux," CoRR, vol. abs/181 1. 1.01457, 2018.
[31] M. Innes, "Flux: Elegant Machine Learning with Julia," " Journal of Open Source Software 2018.
[32] D. P. Kingma and J. Ba, "Adam: A Method for Stochastic Optimization," arXiv, 12 2014.
[33] P. K Mogensen and A. N Riseth, "Optim: A mathematical optimization package for Julia," Journal of Open Source Software, vol. 3, p. 615, 4 2018.
[34] H. Flanders, "Differentiation Under the Integral Sign," The American Mathematical Monthly, vol. 80, p. 615, 1973.
[35] R. Goldman, "Curvature formulas for implicit curves and surfaces," in Computer Aided Geo- metric Design, vol. 22, pp. 632-658, Elsevier, 10 2005.
[36] M. Tao, J. Solomon, and A. Butscher, "Near-Isometric Level Set Tracking," Computer Graphics Forum, vol. 35, pp. 65-77, 8 2016.
[37] M. Fujisawa, Y. Mandachi, and K. T. Miura, "Calculation of Velocity on an Implicit Surface by Curvature Invariance," Information and Media Technologies, vol. 8, no. 4. pp. 674-680, 2013.
[38] M. Kratochvíl, O. Hunewald, L. Heirendt, V. Verissimo. J. Vondrášek, V. P. Satagopam, R. Schneider, C. Trefois, and M. Ollert, "GigaSOM.jl: High-performance clustering and visualization of huge cytometry datasets." GigaScience, vol. 9, pp. 1-8, 11 2020.
11


Reviews:
[{"review_id": "lEfQN--C_H1", "reviewer": "Reviewer_L17k", "summary": "The paper introduces a method for fitting parameterised differential equations to qualitative behaviour: the location and number of bifurcations. Several new techniques are introduced: a loss for matching bifurcations to target locations, a loss for producing the correct number of bifurcations, and finally an approach to differentiating these terms via implicit differentiation methods.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Summary\n\nThis paper is excellent.\n\nThe problem considered is difficult, interesting, and has a clear audience.\n\nThe technical contributions of the paper include several techniques, which to the best of my knowledge are new, technically sound, and frankly quite ingenious.\n\nThe clarity of the presentation is generally good. I was able to understand the paper without too much difficulty. I do have some suggestions for improving this below.\n\nFor these reasons I definitely recommend acceptance.\n\n### Discussion\n\nThe experiments are undoubtedly the weakest part of the paper. Only two are considered (three including the complexity study), and these are on relatively simple synthetic problems. Extensions to less toy, practical, problems would substantially strengthen the paper. In other words, what are some impressive use cases that help sell this technique?\n\nThe differential equations considered are essentially straightforward traditional parameterised differential equations. Do the authors envisage applications to neural differential equations, no doubt on more complicated problems? Perhaps combined with a traditional supervised loss to fit the quantitative behaviour. For example, much of my own work involves fitting neural differential equations to complicated dynamical systems, for which no mechanistic description is available -- only broad qualitative descriptions. It would certainly pique my interest to see an example of this type.\n\nGiven that there is a little extra space at the end of the paper: I think some (brief, intuitive) discussion on the steady-state solver(s) would be the most helpful piece of additional background material to include.\n\n### Minor points\n\nLine 86, equation (5), and elsewhere: $\\frac{d}{ds}|\\frac{\\partial F_\\theta}{\\partial u}|$ does not feature an $s$ to differentiate with respect to. Moreover the definition of the curve parameterised by $s$ (presumably a level set of $F_\\theta$) is not made clear. Leaving these details implicit hurts readability.\n\nEquation (2) is lacking a dividing line between the $p$ and the $\\exists$.\n\nThe upper graphs in Figure 1 could do with a y-axis label $u$.\n\nSection 3.1: it would be helpful to explicitly reference the model (from earlier) which is being optimised. The optimiser used is not stated.\n\nI do not think I would describe the approach as semi-supervised, or that $(|P| - |D|)\\log\\Psi$ represents an unsupervised loss. Whilst these terms are often blurred and without precise meaning: I would take semi-supervised learning to mean that some of the training data could have had a label, but is missing. (For example, training to perform MNIST image classification, except only some 10% of the data has its label.) Likewise I would describe $(|P| - |D|)\\log\\Psi$  as a supervised loss: it pushes together $|P|$ and $|D|$.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No discussion is included, but I think this is fine for a work of this type.", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "bLejGnUOpNw", "reviewer": "Reviewer_cqjt", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well organized and written in a clear and understandable manner. The task that the authors aim to solve is, in general, relevant. Their solution somewhat novel. Previous work on solving this task exists and the authors acknowledge it but don't consider it for direct comparison. Previous work on biasing a goodness of fit with domain knowledge drive the optimization exists, especially coming from the domain of modeling dynamical systems for systems and synthetic biology, however it is not reported in this paper.\n\nThe definition of the task being solved as being semi-supervised might not be appropriate, as it doesn't fit the frame defined by the machine learning community, from where the term is borrowed. The inference of parameters of dynamical systems, given measurements, qualitative or quantitative descriptions of desired or expected behavior better fits a supervised learning task. The task definition and the theoretical analysis are sufficient and convincing.\n\nThe aimed application of the presented work is to inverse problems from biology or engineering. However, the presented applications are toy models and a small synthetic biology model. Inverse problems in biology in general are much more complex than the examples presented in this work and rarely fall in the scope of dynamical systems with codimension one bifurcation. Furthermore, inverse problems from the domain of biology have consistently been proven to be difficult to solve due to being underdetermined, lack of data or various sources of noise. Although the authors acknowledge the model selection problem, they address it superficially in the conclusion section. Additionally, although their minimal models are extensible in the number of parameters they do not analyse their impact on the computational cost. The application of the presented approach to such real world problems is therefore severely limited.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The presented work doesn't have negative societal impact.\n\nThe authors acknowledge the limitations of their work adequately. The main limitation that is not addressed however is the limited applicability of the approach. See main review.", "review_text": "", "extra_scores": {"ethical_concerns": "There are no ethical issues with this paper.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "K735zQLusJ", "reviewer": "Reviewer_DHJo", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "1. The results of this work seem very interesting and important\n2. The paper is well written and clear; but some math parts still need more clarification\n3. I should point out some comments that need to be addressed (please see below):\n\n- In line 79, authors have been mentioned that \"a co-dimension\none bifurcation can be defined by a set of conditions on the determinant of the Jacobian |∂F_θ\\∂u|\". But, based on Sotomayor's Theorem (for codimension-1 local bifurcations), it is only one sufficient condition and to investigate such bifurcations we also need some conditions on \"∂^2F_θ /∂u^2\" and \"∂F_θ\\∂p\" (although they have discussed such conditions briefly in Appendix). So, this sentence is misleading & confusing. \n\n- In Appendix B, it is mentioned \"If we are not specifically looking for transcritical or pitchfork bifurcations, it is sufficient to consider a non-zero Hessian determinant and non-zero vector slope (Eq. (21))\".  But, it is not correct; please see Sotomayor's Theorem for more clarity. Especially, since for saddle-node bifurcations the conditions in Eq.(21) are not enough and such conditions depend on **left eigenvectors** as well. \n\n- Although the authors have mentioned different kinds of bifurcations, from Eq.(2) the type of bifurcation is not clear and they have only discussed relationships between conditions (2) and the conditions for saddle-node bifurcation. Is it possible to add any other conditions from which the type of bifurcation is determined as well? Otherwise, the authors should clarify that (2) is related to only one bifurcation e.g. saddle-node bifurcation.\n\n- In line 91, it is mentioned \"Pitchfork bifurcations are special cases of the saddle-node where a single steady state splits into two stable and one unstable steady state (Figure 1B).\" But, this is mathematically **incorrect** since pitchfork bifurcations are **not** special cases of the saddle-node bifurcations; please see Sotomayor's Theorem for more clarity. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers agree that the paper propose an interesting approach to the problem of estimating parameters in differential equations. Although some reviewers have some technical concerns at their first reviews, they satisfy the authors' responses to address those. It seems that, although there are some points that should be modified from the current form including the terms used in the paper (such as semi-supervised learning), I think we can expect the authors modify the paper in the camera-ready by reflecting the discussion. Based on these, I recommend acceptance (poster) for this paper.
Decision: Accept (Poster)