Geometric analysis of mixed-mode oscillations in a model of electrical activity in human beta-cells

Human pancreatic beta-cells may exhibit complex mixed-mode oscillatory electrical activity, which underlies insulin secretion. A recent biophysical model of human beta-cell electrophysiology can simulate such bursting behavior, but a mathematical understanding of the model’s dynamics is still lacking. Here we exploit time-scale separation to simplify the original model to a simpler three-dimensional model that retains the behavior of the original model and allows us to apply geometric singular perturbation theory to investigate the origin of mixed-mode oscillations. Changing a parameter modeling the maximal conductance of a potassium current, we ﬁnd that the reduced model possesses a singular Hopf bifurcation that results in small-amplitude oscillations, which go through a period-doubling sequence and chaos until the birth of a large-scale return mechanism and bursting dynamics. The theory of folded node singularities provide insight into the bursting dynamics further away from the singular Hopf bifurcation and the eventual transition to simple spiking activity. Numeri-cal simulations conﬁrm that the insight obtained from the analysis of the reduced model can be lifted back to the original model.


Introduction
Glucose-induced insulin secretion from pancreatic betacells involves electrical activity, which leads to opening of voltage-gated Ca 2+ channels and exocytosis of insulin-containing granules [1,12,24]. Mouse beta-cells provide the typical experimental material for the study of electrical behavior underlying insulin release, and decades of studies have shown that the membrane voltage of mouse beta-cells exhibit bursting oscillations composed of active phases of action potential firing from a plateau interspersed by silent, hyperpolarized phases. This behavior, and complex variations of it, are believed to underlie pulsatile insulin release [4,12,28]. Mathematical modeling has played an important role in unraveling the mechanisms underlying the complex electrophysiological behavior of mouse beta-cells [4,19,25], and these models represent a prototypical example of so-called square-wave bursters that can be studied as a slow-fast system using time-scale separation [3,22].
More recently an increasing number of studies have been performed on human beta-cells, which show similarities but also important electrophysiological differences to their counterpart from mice [5,16,21,24]. Human beta-cells typically show spiking or very fast bursting, with about half of the cell population exhibiting each type of behavior. In contrast to square-wave bursting observed in mouse beta-cells, bursting in human beta-cells is often characterized by an active phase with very few and small oscillations around a depolarized plateau [2,5,16,23]. This behavior resembles activity in certain pituitary cells and has been coined pseudo-plateau bursting [26].
A decade ago the first mathematical model of human beta-cell electrophysiology was developed [20]. This model and its descendants [14,15,17,23] can simu-late both spiking and fast bursting activity. However, so far no mathematical analysis of this model, in particular with respect to the dynamics underlying pseudoplateau bursting, has been performed. Inspired by analysis of models exhibiting pseudo-plateau bursting in pituitary cells [26,29], we speculate that mixed-mode oscillations (MMOs) due to folded singularities orchestrate pseudo-plateau bursting in models of human betacells. We refer to MMOs as complex patterns that arise in dynamical systems, in which Small Amplitude Oscillations (SAO) and Large Amplitude Oscillations (LAO) are interspersed [3,6,8].

Theoretical background
In the following sections we will reduce the original model of electrical activity in human beta-cells [20] to a three-dimensional model that will allow us to apply MMO theory. In this section we outline this theory, mainly following Desroches et al. [8].
We consider the three-dimensional fast-slow system of singularly perturbed ordinary differential equations ẋ = f (x, y, z) , y = g 1 (x, y, z) , z = g 2 (x, y, z) , (1) with f, g 1 , g 2 sufficiently smooth functions and a small parameter 0 < 1. Let t denote the independent variable of (1) andẋ the derivative of x with respect to t. In this system x is the fast variable and y, z are the slow variables, and t is referred to as the slow time scale.
By switching to the fast time scale τ = t/ , one obtains the equivalent system where x denotes the derivative of x with respect to τ . As → 0 the trajectories of (2) converge during fast epochs to solutions of the fast subsystem or layer equations implying that the fast flow is almost horizontal, following the principal direction of the fast variable x.
Conversely during slow epochs trajectories of (1) converge to solutions of the slow subsystem or reduced problem: which is a system of differential-algebraic equations.
Geometric singular perturbation theory (GSPT) is used to understand the dynamics of the full systems (1) and (2) for > 0 from the (singular) fast and slow subsystems (3) and (4). The algebraic equation in (4) defines the so-called critical manifold and the slow flow determined by (4) is constrained to lay on this manifold. S is supposed to be a smooth manifold. Note that the points of S are equilibrium points for the layer equations (3). GSPT guarantees that normally hyperbolic invariant manifold S persists as locally invariant slow manifolds S of the full problem (1) for sufficiently small . Moreover the restriction of the flow of (1) to S is a small smooth perturbation of the flow of the reduced problem (4) [9]. Normal hyperbolicity is often difficult to verify when there is only a single time scale. However, in our slow-fast setting, S consists entirely of equilibria (of the layer equations (3)) and therefore the requirement of normal hyperbolicity of a compact subset S 0 ⊂ S is satisfied as soon as all p ∈ S 0 are hyperbolic equilibria, i.e., if f x := ∂f /∂x is uniformly bounded away from the imaginary axis for all p ∈ S 0 . In that case the critical manifold S 0 is normally hyperbolic.
GSPT breaks down at points on the critical manifold where normal hyperbolicity fails. This happens at points on S where their projection onto the space of slow variables is singular, i.e., on the fold lines L, defined as the set of equilibria of (3) with a zero eigenvalue, The stability of points p ∈ S as equilibria of the layer problem (3) depends on the sign of f x : the fold lines divide the critical manifold S into attracting sheets S a (where f x < 0) and repelling sheets S r (where f x > 0). The critical manifold S is assumed to be (locally) a folded surface, i.e., and we suppose f xx = 0 for every point on the fold line L, i.e., This condition ensures that S has a fold, and not a deflection tangent, in L. Without loss of generality, we assume that f z (p)| p∈L = 0 holds locally. It follows by the implicit function theorem that S is locally given as a graph z = φ(x, y). Hence the reduced problem (4) can be projected onto the (x, y)-plane, yielding where the first equation is obtained by differentiating the algebraic equation f = 0 (which generates the folded surface) with respect to t. Indeed, this equation is no longer necessary after imposing the equivalent constraint z = φ(x, y); conversely its derivative allows us to make explicit the evolution of the fast variable on the slow manifold S. Finally, desingularization (i.e. rescaling time by dt = −f x dt ) removes the singular term on the fold and gives the desingularized systeṁ . (9) The desingularized system (9) is equivalent to the original one (8) on the attracting sheets S a but has opposite orientation on the repelling sheets S r due to the time rescaling: for every point p ∈ S a we have by definition f x (p) < 0, hence rescaling time by the factor −f x (p) > 0 does not change the orientation; conversely, −f x (p) is negative on S r and thus changes the direction of the flow.
Singularities of the desingularized system (9), i.e., points whereẋ =ẏ = 0, can be classed as ordinary or folded. Ordinary singularities are true equilibria on S of the initial system (1) where g 1 is vanishing (and g 2 = 0 as well, since f y and f z are nonvanishing in the vicinity of L). Folded singularities, on the other hand, are points along the fold lines L (hence where f x = 0) where f y g 1 + f z g 2 | z=φ(x,y) equals zero. These are characterized as folded saddles, nodes or foci based on the eigenvalues σ s and σ w of the Jacobian J of (9) evaluated at the folded singularity. We use notation so that |σ s | ≥ |σ w | and denote σ s and σ w the strong and weak eigenvalue, respectively.

Canards
The trajectories of the slow flow that lie along the eigendirections of a folded saddle, or a folded node, connect the two sheets of the critical manifold through the folded singularity in finite time. In these points the slow flow switches from incoming to outgoing. Such solutions are called singular canards and their persistence under small perturbations can give rise to complex dynamics. We call the singular canard corresponding to the strong (respectively weak) eigendirection the strong (respectively weak) singular canardγ s (respectivelyγ w ).
For the case of a folded node, the strong singular canardγ s and the fold curve L bound sector of trajectories that cross from S a to S r by passing through the folded node. This sector is called the funnel of the folded node, and the funnel contains necessarily the weak singular canard. In the system we are going to study, the relevant folded singularities are attracting folded nodes.
Singular canards persist under small perturbations > 0, and their perturbation are simply denoted canards. In geometric terms, a canard corresponds to the intersection of the manifold S a, and S r, extended by the flow to the vicinity of the folded singularity [8]. Canards are not unique since the corresponding invariant manifolds S a, and S r, are not unique. For a fixed choice of the invariant manifolds we call their intersections maximal canards. Maximal canards can be seen as perturbations of singular canards that preserve the original behaviour. The singular canardγ s (the strong canard ) always perturbs to a maximal canard γ s , and in general, the singular canardγ w (the weak canard ) also perturbs to a maximal canard γ w [6,8]. We call γ s and γ w primary canards. In addition to these, there are other maximal canards, which we call secondary canards. The presence of secondary canards divides the sector between primary maximal canardsγ w andγ s into subsectors. Furthermore, trajectories within each subsector will be forced to take a certain number of turns aroundγ w when passing through the folded node resulting in SAOs [6,8].

MMOs in a model of human beta-cells
Our starting point is a slightly modified version of the Hodgkin-Huxley type model of electrical activity in human beta-cells [20]. Briefly, this seven-dimensional nonlinear model represents the dynamics of the beta-cell membrane potential V , which is controlled by currents (normalized to cell capacitance) through a set of ion channels. Channel gating operating on different time scales makes the dynamics highly nontrivial and different kinds of dynamical patterns can be observed depending on the model parameters. In the following we refer to this model as the 7d-model. Further details and biophysical motivation can be found in the original paper [20]. For completeness we provide the equations of with currents and steady-state activation (m X,∞ ) and inactivation (h X,∞ ) functions of the form Finally, the voltage-dependent time scale of m K V is The model reproduces the three main types of electrical behavior seen in human beta-cells: silence characterized by the absence of action potential generation, spiking activity consisting of individual action potentials, and pseudo-plateau bursting where small action potentials fire from a depolarized plateau (Fig. 1). Table 1 Characterization of variables according to their time scales. Note that τ mKV (V ) depends on V , and the range of its values is given.
Changing the maximal conductance of the hERG potassium current, g hERG , suffices to switch between the different types of behavior. In the following we are mainly interested in understanding the dynamics of pseudoplateau bursting.

Model reduction
To apply the theory exposed in Section 2 we reduce the initial 7d-model to a 3-dimensional system. Based on the time scales shown in Table 1 we divide the variables into three categories (fast, medium and slow). The (fast) variable V is excluded from this analysis since all the other variables depend on it.
For the fastest variables we use quasi steady-state approximations as frequently done for Hodgkin-Huxley types models. From a mathematical point of view, this implies that these variables should be approximated by their steady-state functions. Simulations shown in Fig. 2 confirm that the quasi steady-state approximation is reasonable for the three fast variables h CaT , h N a , m BK .
The medium-fast variables h CaL and m KV are not well approximated by their steady-state functions and we use another widely applied approach due to Fitzhugh [10]. Since the two variables have similar timescales, we plot m KV versus h CaL , and fit the curve with a linear function. As seen in Fig. 3, for the case of bursting, which is our main interest, an adequate fit for the active phase is given by We thus obtain the three-dimensional system (from now on we use h and m in place of h CaL and m hERG , respectively), We denote this system the 3d model. Simulations show that the three types of dynamics (hyperpolarization, spiking, bursting) present in the original 7d model are preserved in this model (Fig. 4).

Geometric singular perturbation theory
The characteristic time constant for the V equation in (15) is the reciprocal of the maximal conductance normalized to cell capacitance, i.e., τ V = 1/ max{g X } = 1/g N a = 2.5 ms.
We conclude that V is the fastest variable and, after scaling of time, rewrite the system in the form studied in Section 2, witht = t/τ m , = 1/(g N a τ m ) = 0.025 andĨ = I/g N a , This allows us to apply GSPT to our 3d-model to understand the mechanism organizing pseudo-plateau bursting, which is observed for example with g hERG = 0.16 nS/pF (Fig. 4). We remark that g 2 (respectively g 1 ) does not depend on h (respectively m) and that the dependence on m  (respectively h) is linear. In addition, f is linear with respect to both h and m, and has the form for appropriate functions A, B, C. Note that the form (18) is rather general for Hodgkin-Huxley-type, conductancebased models. It follows from the above that when we solve f = 0 in order to calculate the folded surface S defined in (5), we can find m as an explicit function of V and h, i.e., Fold curves are obtained by asking f V = 0, which yields after imposing (19). Note that the fold curve does not depend on g hERG .
We proceed with the analysis by studying system (17) in the limit case when → 0. We project the system onto the critical manifold, described by the (V, h) coordinates, and find the slow desingularized problem of the form of (9), evaluated at m = φ(V, h). Using (17), (18) and (19) this system becomeṡ For standard parameters and g hERG = 0.16 nS/pF, the desingularized system (22) has three equilibria at points whereV = 0 and at least one of f V and g 1 is vanishing. Two of these are plotted in Fig. 5, while the third one is out of the range of h ∈ (0, 1); this latter is a folded focus, and hence it is of no interest for the analysis of canards [8] and for the dynamics of the system.
The other two equilibria of the desingularized system are an attractive folded node and a saddle. The saddle is a real equilibrium of system (17) since it is not on the fold and g 1 = 0 (and g 2 = 0 as well, see Section 2). On the other hand, the attractive folded node lies on one of the fold curves and therefore represents the point where the canard trajectories go from the attractive manifold S a to the repelling one S r (see Fig. 5). The singular canards and the consequent funnel (shadowed gray region) limited by the strong singular canard and the fold curve on the attractive manifold S a are shown in Fig. 5, together with the simulation previously computed in Fig. 4, projected onto the 2d-plane (V, h CaL ). The trajectory enters the funnel and starts oscillating around the weak singular canard passing through the folded node. The system is attracted towards the saddle equilibrium and is then, after a few oscillations, pushed away from the saddle in the unstable eigendirection. Finally the global return mechanism yields a LAO and permits the system to repeat the cycle.

Maximal Canards
We now go into greater details for the study of the system for > 0. Recall that S a and S r are perturbed to attractive and repelling slow manifolds S a, and S r, , and that maximal canards appear in the vicinity of the folded node as the intersections of S a, and S r, . Further, secondary canards split the funnel -more precisely they split the sector of secondary canards -into subsectors. When the solution of the system enters the funnel, it enters in one of these specific subsectors, which are determined by the primary and secondary maximal canards.  Fig. 5 Analysis of the desingularized system and simulations of the 3d model using default parameters, and g hERG = 0.16 nS/pF. Fold curves are plotted in green, and the attracting and repelling sheets are indicated by S a and S r , respectively. The blue point indicates the folded node while the red one is the ordinary saddle equilibrium. The eigenvector directions of the two singularities are plotted with the corresponding colors, and using dashed thick and dotted lines for the strong and weak eigenvectors, respectively. In particular, the strong (γ s ) and weak (γ w ) singular canards correspond to the blue dashed thick and dotted lines, respectively. The funnel of the folded node is indicated by the shaded grey region. The projection of the simulation of the 3d model is shown by the black curve.
We find S a, (respectively S r, ) by evolving the 3dmodel forward (respectively backward) in time starting from points of a line on S a (respectively S r ) at some distance from the fold. Such points are -close to the slow manifold, and are rapidly attracted to S a, (respectively S r, ) in forward (respectively backward) time. Our approach can be seen as a simplified alternative to finding the intersections of the slow manifolds by defining an appropriate Boundary Value Problem [7]. Figure 6A shows these geometric structures for g hERG = 0.16 nS/pF, corresponding to the dynamics simulated in Fig. 4. As predicted by the theory, S a, and S r, are twisted and intersect in twisting curves (maximal canards). This fact is more clearly seen when plotting the intersection between S a, , respectively S r, , and a plane defined by m being constant (Fig. 6C). These intersections define two curves in the plane, and their crossings correspond to the maximal canards. The solution to the 3d-model (black) enters the region near the folded node along the attracting sheet in the section between the first (γ 1 ) and the second (γ 2 ) secondary canards, and hence perform two SAOs before leaving the fold region. The return mechanism, which yields a large amplitude oscillation, then reinjects the system into the funnel and the dynamics is repeated, resulting in 1 2 mixed-mode oscillations. In contrast, for g hERG = 0.18 nS/pF, the return mechanism injects the system into the section between the strong canard γ s and first secondary canard γ 1 (Fig. 6BD). Therefore the system performs only one SAO during the passage through the fold region and thus exhibits 1 1 MMOs.

Bifurcations
We now aim to understand how changes in g hERG lead to different dynamical patterns. The relation between g hERG and the folded singularity can be found from imposingV = 0 in (22) with h(V ) = Ψ (V ) and solving for g hERG . Similarly, the curve of the full system equilibrium is found by isolating g hERG from dV dt = 0 in (17) with m = m ∞ (V ) and h = h ∞ (V ), i.e., by solving φ(V, h ∞ (V )) = m ∞ (V ). As shown in Fig. 7 these curve intersect at approximately (0.06 nS/pF, −25.3 mV), and thus a folded saddle-node (FSN) of type II [8] is present at these values. For all other values of g hERG , the folded singularity is a saddle. For g hERG < 0.06 nS/pF, the full system equilibrium is stable, corresponding to a permanently depolarized state, and located on the attracting sheet of the slow manifold.
It is known that a singular Hopf Bifurcation is often occurring near FSN of type II [8,11]. Indeed, by calculating the eigenvalues of the full-system equilibriums, we find a Hopf bifurcation (HB) at g hERG = 0.1095 nS/pF. This HB is supercritical. We performed a scan over a range of g hERG values to follow the limit cycle born in the HB. Fig. 7 shows the numerically calculated maximum and minimum values. We see that the system goes through a period-doubling cascade as g hERG increases, entering chaos, which is interrupted at g hERG ≈ 0.1537 nS/pF where a sharp transition of the lowest minimum occurs. This bifurcation, which leads to MMOs, is most likely due to the intersection of the unstable manifold of the saddle-focus and the repelling sheet of the slow manifold [11,18].
As g hERG increases further, the MMOs become regular 1 2 oscillations, followed by a transition to 1 1 MMOs, and finally regular 1 0 spiking with no SAOs, as indicated in Fig. 7. These transitions can be understood from the return mechanism re-injecting the system into different sectors of the folded-node funnel (see Section 3.2.1 and Fig. 6 regarding g hERG = 0.16 nS/pF versus g hERG = 0.18 nS/pF) until it misses the funnel and relaxation-type spiking oscillations appear. Thus, as g hERG increases, the SAO-generating mechanism shifts from being based on dynamics associated with the singular HB to being organized by canards and the twisted slow manifolds around the folded node. In addition, the maximal number of canard-induced SAOs, which can be calculated from the ratio µ of the eigenvalues of the folded node as [6,8] where · denotes the floor function, decreases as g hERG increases (Fig. 7, lower panel).

Discussion
Bursting electrical activity has been proposed to evoke more hormone release than simple action potential firing [27], and it is therefore important to understand the mechanisms that drive this kind of electrical behaviour. Pseudo-plateau bursting, as observed in human betacells [2,5,16,23] and reproduced by the mathematical model studied here [20], can be seen as MMO dynamics and thus studied with the theory of folded singularities as done in this work.
To perform the mathematical analysis we performed a series of approximations to reduce the system to a 3d-model, which exhibited the same kind of qualitative patterns as the 7d-model. It is natural to ask whether the analysis performed for the 3d-model can be lifted back to the original model. Indeed, as shown in Fig. 8, the simulated dynamics of the 7d-model appears to be organized by the structure found for the 3d-model. When projected onto the (V, h Ca , m hERG ) space, the 7d-system approaches the folded node found for the 3dmodel and performs SAO around the weak canard. It then approaches the saddle point located on the repelling sheet before spiralling away. A global return mechanism takes the orbit back to the attracting sheet and the cycle is repeated.
In simulations of the original model [20], as well as in the reduced 3d-model presented here, qualitative changes from spiking to bursting to silent electrical behavior are seen as the hERG conductance g hERG is reduced ( Figs. 1 and 4). The same kind of qualitative changes can be seen for other modifications to the parameters of the potassium currents in the model [20]. Based on the analysis performed here, we propose that these changes can be understood as follows. Starting from a permanently depolarized state, increasing (hyperpolarizing) potassium current moves the equilibrium from the attracting to the repelling sheet of the slow manifold, passing through a FSN of type II. As the potassium conductance, here g hERG , increases further, the depolarized equilibrium loses stability in a singular Hopf bifurcation. Following intervals with regular and chaotic oscillations, bursting in the form of MMOs appear when the chaotic attractor is destroyed when the

Bifurcations
We now aim to understand how changes in g hERG lead to di↵erent dynamical patterns. The relation between g hERG and the folded singularity can be found from imposingV = 0 in (22) with h(V ) = (V ) and solving for g hERG . Similarly, the curve of the full system equilibrium is found by isolating g hERG from dV dt = 0 in (17) with m = m 1 (V ) and h = h 1 (V ), i.e., by solving (V, h 1 (V )) = m 1 (V ). As shown in Fig. 7 these curve intersect at approximately (0.06 nS/pF, 25.3 mV), and thus a folded saddle-node (FSN) of type II [7] is present at these values. For all other values of g hERG , the folded singularity is a saddle. For g hERG < 0.06 nS/pF, the full system equilibrium is stable, corresponding to a permanently depolarized state, and located on the attracting sheet of the slow manifold.
It is known that a singular Hopf Bifurcation is often occurring near FSN unstable manifold of the saddle-focus and the repelling sheet of the slow manifold intersect. As the potassium conductance increases further the MMOs are organized mainly by the folded node, and when the return mechanism no longer inject the orbit into the funnel of the folded node, relaxation-type spiking oscillations appear. Eventually, the equilibrium moves to the other attracting sheet of the slow manifold and the system becomes stable and hyperpolarized.
Although our analysis of the 3d model appears to capture the dynamics of the 7d model satisfactorily, some discrepancies may be observed. For example, the 7d model may exhibit 1 3 MMOs (Fig. 1), which were not observed in the 3d model when scanning through different g hERG values (Fig. 7). This is likely due to the model simplifications. Indeed, the approximation (14) is not very accurate, and misses for example the final oscillation of the active phase at (h CaL , m KV ) ≈ (0.4, 0.14), see Fig. 3. Future studies could analyze the 4d model  Fig. 7 Overview of the behavior of the 3d-model when g hERG changes. For low values of g hERG the model has stable equilibriums indicated by the full curve. At g hERG ≈ 0.11 nS/pF the system undergoes a supercritical Hopf bifurcation (HB) and the equilibrium becomes an unstable saddle-focus (dotted curve). The dash-dotted line indicates the location of the folded node. This curve intersects the branch of equilibriums at g hERG ≈ 0.06 nS/pF where a folded saddle-node (FSN) of type II is present, as shown in the upper inset (corresponding to the gray horizontal strip in the main figure). The blue points show maxima and minima for the attractive orbits originating from the HB. The stable limit cycle born in the HB undergoes a period-doubling (PD) cascade to chaos as shown in the lower inset (corresponding to the gray vertical strip in the main figure). At g hERG ≈ 0.1537 nS/pF, an abrupt transition to MMOs occurs, as indicated by the gap in the diagram in the lower inset (red triangle). Following transitions through windows of apparent chaotic behavior, regular 1 2 and 1 1 MMOs are present for large intervals of g hERG , as indicated, until only 1 0 large-scale oscillations corresponding to spiking appear for g hERG > 0.21 nS/pF. The lower panel shows the maximum number of SAOs calculated from (23) as a function of g hERG .
(V, h CaL , m KV , m hERG ) as a 1 fast -3 slow variable system. Concerning the 3d model, it may be seen as a three timescale problem [13], in particular in case the time-constant for m hERG is increased compared to the (biophysically realistic) value used here. Such an analysis is left for future work. Concluding, we have shown an example of how advanced mathematical analysis based on geometric sin- Fig. 8 The analysis for the 3d-model provides insight into the 7d-model. A simulated orbit (black) of the 7d-model (from Fig. 1) is projected onto (V, h Ca , m hERG ) space and compared to the geometric structures found for the analysis of the 3d-model with g hERG = 0.16 nS/pF: The critical manifold composed of the attracting (S a ) and repelling (S r ) sheets separated by the fold curves (green curve) is shown with the folded node (blue point) and the saddle (red). The computed orbit shows small oscillations starting near the folded node that spiral around the weak singular canard (blue line) before escaping towards the other fold curve. A relaxation-like return mechanism takes the system back to the attracting sheet S a , yielding the first peak seen in the simulation in Fig. 1. S a is then followed towards the upper fold and the events repeat. gular perturbation theory can provide insight into a complex but biologically realistic mathematical models of electrical activity underlying insulin release in man.

Funding
This work was supported by MIUR (Italian Minister for Education and Research) under the initiative Departments of Excellence (Law 232/2016). The funding source had no role in study design; in the collection, analysis, and interpretation of data; in the writing of the report; or in the decision to submit the paper for publication.

Conflict of interest
The authors declare that they have no conflict of interest.

Code availability
The computer code used in this work is available at http://www.dei.unipd.it/~pedersen/ or from the corresponding author upon reasonable request.

Author contributions
M.G. Pedersen developed the study conception and design. Both authors performed mathematical analysis, developed computer code and prepared figures. The draft of the manuscript was written by M.G. Pedersen, and S. Battaglin commented on previous versions of the manuscript. Both authors read and approved the final manuscript.