Fractional critical slowing down in some biological models

The critical slowing down (CSD) phenomenon of the switching time in response to perturbation β (0 < β < 1) of the control parameters at the critical points of the steady state bistable curves, associated with two biological models (the spruce budworm outbreak model and the Thomas reaction model for enzyme membrane) is investigated within fractional derivative forms of order α (0 < α < 1) that allows for memory mechanism. We use two definitions of fractional derivative, namely, Caputo’s and Caputo-Fabrizio’s fractional derivatives. Both definitions of fractional derivative yield the same qualitative results. The interplay of the two parameters α (as memory index) and β shows that the time delay τ D can be reduced or increased, compared with the ordinary derivative case (α = 1). Further, τ D fits: (i) as function of β the scaling inverse square root formula 1 / β at fixed fractional derivative index (α < 1) and, (ii) as a function of α (0 < α < 1) an exponentially increasing form at fixed perturbation parameter β.


Introduction
Bistable systems in many branches of sciences (physics, biology, . . . ) and engineering are characterized by the co-existence of two stable states, where the system switches from one stable state to other state by means of changing one or some of the system control parameters [1][2][3][4]. The associated transient phenomena of lengthening the switching time between these two stable states, called critical slowing down (CSD), happens upon perturbing one of the parameters at the critical (switching-on or -off) points of the charactertistic bistable curve [5][6][7][8]. It has been suggested that, CSD may serve as a universal indicator of how a complex physical system (such as brain, ecosystems, climate and financial markets) approaches a threshold [9][10][11][12], and as well serving as an indicator of transitions in two-species biological models, which exhibit Hopf bifurcation or hysteresis transition [13]. For our specific current concern, the CSD phenomenon has recently been investigated by us in [14] for some biological bistable models, namely.
Specifically, our investigation in [14] was concerned with the nature of transition between the two stable states, and the verification of the inverse square root scaling law, for the switching time delay (τ D ) at the critical switching-on and -off points, independent of the type of nonlinearity in the model rate equations. The model rate equation in model a) is of first order ordinary differential equation (ODE), while in model b) the model rate equations are coupled first order ODEs.On the other hand, fractional calculus, a field of mathematics that deals with the analysis of derivatives and integrals of fractional (or even complex) order, has its applications in diverse areas of science and engineering. The associated fractional differential equations (FDEs) are widely and successfully used in mathematical modelling in a variety of fields. We refer the reader to the extensive list of major works and applications in the area of fractional calculus cited in ( [17][18][19][20] and refs. therein). In ordinary calculus, the first order derivative of a function f(t), namely f′(t) df dt is the instantaneous rate of change of f(t) over the infinitesimal time period, t → 0, that is, local time effect. In fractional calculus, the physical meaning is non-local, as the time domain is manifested as a memory (or time delay) effect and the current state of the system depends on its earlier history. Moreover, in fitting with test data of various models of memory phenomena, the order of the fractional derivative serves as an index of memory [21,22]. FDEs of arbitrary real order are not in general easy to solve analytically [23]. However, the numerical method based on Laplace transform technique is a basic one and applicable for a wide class of initial value problems for FDEs, [23][24][25][26]. Recent fundamental computational methods are found in [27,28]; and refs, therein.
Experimentally speaking, fractional derivative models (FDMs) are in excellent agreement with experimental data in many branches of science and engineering. Two specific examples we quote. Further, CSD or more generally instability mechanism and chaos, have been investigated at large in fractional order dynamical systems in fields, like, fluid flow [31-35], neurology and biological phenomena ( [36][37][38]; refs. therein) to account successfully for memory (timedelay) and special non-local effects. For example.
1. The Landau model that describes the fluid flow from laminar to turbulent has been examined within a fractional rate equation model [35] in order to account for memory effect. This transition to turbulence due to CSD shows that the turbulent fluctuations depend on memory of inverse power law decay in agreement with experiment [39]-slower than in the case of no memory (ordinary derivative case) of turbulent fluctuations decaying exponentially, 2. Capacitive memory due to fractional order cardiomyocyte dynamical model [37] alters the electrical signaling in cardiac cells in a manner that promote or suppress electrical instability (known as alternans). 3. The use of a fractional order mathematical model to study the signaling process in nerve cells (like, neuron) due to incorporated strong memory effects [36] has been interpreted as a neuronal disorder (Parkinson disease).
The concern of the present paper is to adopt the corresponding FDEs in both models a) [3,4,15] and b) [4,16], referred to above, in order to incorporate for memory effects and examine effects of the fractional derivative order parameter (α), (0 < α < 1) on the time delay (τ D ) associated with the CSD phenomena examined in the no-memory case [14]. We use and compare two definitions of fractional derivatives, namely, Caputo's [40] and Caputo-Fabrizio's [21,22] definitions. Both definitions have the advantage of dealing with initial conditions of the variables and their integer derivatives suitable in most physical processes, like models a) [3,4,15] and b) [4,16] referred to above. As a main result, it is found that Caputo's and Caputo-Fabrizio's definitions of fractional derivatives yield the same qualitative results of reduced time delay τ D at fixed perturbation of the concerned control parameter, with smaller values of the fractional derivative order α. The small quantitative difference in τ D is due to the different convoluted kernels (that model the memory or delay effect) in [21,22,40].This paper is presented as follows. In section 2), we present the model differential equations in both ordinary and Caputo's fractional derivative forms, for both models. In section 3), we present the computational results for the transient switching. Section 4) presents a summary of the results. In Supplementary Appendix A, a brief background of the model ODEs (eqa (1) and. 2) below) representing the two biological models referred to above is given, while Supplementary Appendix B presents a guideline for Euler's numerical method to solve fractional FDE.

The model equations
Here, we first present the model DEs of the two biological models (the Spruce-budworm and Thomas reaction models) in their ordinary derivative forms. (A brief background of these model ODEs are given in Supplementary Appendix A). Second, we present the corresponding fractional derivative forms, according to the two formulations of Caputo's [40] and Caputo-Fabrizio's definitions [21].
where N(τ) is the budworm's population, τ = rt is normalised time, r is the linear birth rate and K is the constant carrying capacity which is related to the foliage (food) available on the trees in the absence of birds. The constant F = p o A/r is the predation population with rate p o and A is the (positive) predator attack rate and B is the threshold measure of the budworm population. The predation will approach an upper level value, lim N→∞ FN 2 /(N 2 + B 2 ) F as N increases.

The Thomas reaction model
The mechanism of this model is based on a basic reaction in an enzyme membrane, between the substrate oxygen and uric acid. The model equations of the system in a dimensionless form are [4,16]: ) .

(2b)
Frontiers in Physics frontiersin.org Here, u and v represent the uric acid and the oxygen being supplied at constant rates a and γb, respectively, where, a, ℓ, k, γ and b are all positive real constants. The factor u(τ)v(τ)/ (1 + u(τ) + ku 2 (τ)) exhibits substrate inhibition: it increases (decreases) when u is small (large), with measure of inhibition's severity equal to k.
In [14], the model Equations 1, 2 were analysed in detail (theoretically and computationally) regarding regions of bistability, the CSD phenomenon at the critical (switch-up and -down) points of the bistable curves and the verification of the inverse square root scaling law of the switching time delay [7,41].

Fractional derivative cases
In this case, Equations 1, 2 take the following forms; and, respectively, where d α dτ α denotes the fractional derivative of order α (0 < α < 1). There is no unique definition of fractional calculus (FC), derivatives and integrals. Definitions of FC are too many and still -up to date -increasing. Here, we use and compare two definitions of fractional derivatives of a continuous function f(τ) on (0, τ), namely, Caputo's [40] and Caputo-Fabrizio's [21] derivatives.

Transient switching and time delay
The switching time at the critical (switch-on and -off) points of the characteristic steady state bistable curves (N vs K) according to the FDE 3), or (u and v vs a) according to the FDEs 4) with both Caputo's and Caputo-Fabrizio's fractional derivatives, Eqs. 5 and 6, respectively, are investigated by solving these FDEs numerically using the fractional Euler's method developed in [28,48] (see Supplementary Appendix B for guidelines). This is done by replacing the control (input) parameter K in Equations 1-3) by K c ± β, or a in Equations 2, 4 by a c ± β, where β (0 < β < 1) is a small real perturbation of the relevant control parameter, and K c , a c are the initial (switch-on or switch-off) points of the bistable curves. Results are compared with the ordinary derivatives case (α = 1) [14].

The spruce budworm model
The switching-on and off -points, A on and A off , respectively, of the steady state bistable curve (N vs. K) according to the ODE, Eq. 1, or the FDE; Eq. 3, i. e., dN dτ d α N dτ α 0, are shown in Figure 1, for fixed values of the parameters F and B (c.f [14]). For fixed positive perturbation parameter β = 0.1, the time delay τ D to switch up to the upper branch of the bistable curve, Figure 2, is reduced in both cases of the fractional derivatives with smaller values of α, (0 < α < 1), compared with the ordinary derivative case (α = 1). This is confirmed in Figure 3 where for fixed 0 < β < 1, τ D vs. α best fits exponentially increasing function for α ∈ (0, 1) in both cases of fractional derivatives. Note in Figure 2, τ D is slightly reduced in the case of Caputo-Fabrizio's fractional derivative, compared with Caputo's fractional derivative case. Further, for fixed fractional parameter α = 0.25, for example, the lesser the perturbation parameter β, the larger is τ D (i.e. slowing down)-Figures 4-like the case of α = 1 [14]. For fixed negative value of perturbation-Figures 5-at the switching-off point A off (in Figure 1), we have the same qualitative behaviour as in Figure 2, but with smooth delayed switching to the lower branch.
In both cases of positive and negative perturbations β) at the switching-on and -off points, A on and A off , respectively, the time delay

FIGURE 5
Data as in Figure 2, but with negative perturbation, K c − β at A off , where K c =3.0199.

The Thomas reaction model
The steady state bistable curves for the Oxygen and uric acid concentrations u, v, respectively, against the supplied rate a, according to Eq. 2 or 4) are shown in Figure 7, for fixed values of other system parameters [14]. For positive perturbation β in the ordinary derivative case (α = 1) at the switching-on point, A on in Figure 7, the transient oxygen concentration u(τ), Figure 8, shows similar qualitative behaviour of reducing τ D in both cases of Caputo's and Caputo-Fabrizio's fractional derivatives, but with smaller quantitive difference. The same behaviour occurs with negative perturbation at the switching-off point A off in Figure 7. Similar qualitative behaviour is also exhibited for the transient uric acid concentration v(t) for α = 1 [14] and α < 1. The time delay τ D in both cases of u(τ) and v(τ) against the fractional parameter α and the perturbation parameter β shows similar qualitative behaviour as in Figures 3, 6, respectively.

FIGURE 7
The steady state bistable curves, u and 0.12v, versus the control parameter, a, for fixed parameters K =20, B =100, γ = l =1. Fractional order mathematical models generalise the concept of ordinary differentiation to incorporate memory (time delay) and spatial non-local effects, and hence provide extra fractional parameters to interpret/predict the dynamical behaviour of the concerned model and capture more of its details.In this paper, we have investigated the switching time response at the critical switching-on and -off points of the bistable curves related to two biological models, namely, the spruce budworm outbreak model [3,4,15] and the Thomas reaction model for enzyme membrane [4,16] within fractional order models. Two definitions of fractional derivatives of order α, (0 < α < 1) have been used, namely, Caputo's [40] and Caputo-Fabrizio's [21,22] fractional derivatives. Our study shows the following.
(i) The two definitions use convolution kernels of different variability that model the memory effect, namely, as power function [40] and as exponential function [21]. Both definitions yield the same qualitative results, (ii)-(iv) below, for the two biological models referred to above. The small quantitative variance in the results is due to the different mathematical forms for the memory or delay effect. (ii) The switching time τ D due to the perturbation in the control (input) parameter, at the critical points of the bistable curves, is reduced further in the fractional derivative case (0 < α < 1), compared with the ordinary derivative case (α = 1) [14], (iii) For fixed perturbation β, τ D as a function of the fractional derivative parameter, α, (0 < α < 1) fits an exponential form, i.e., τ D is reduced with strong memory index (α ≪ 1) and, (iv) The switching time τ D as a function of the perturbation parameter β fits the scaled inverse square root law 1 β √ at fixed fractional derivative index (α < 1) as in the ordinary derivative case (α = 1) [14]. This is a further indication of the universality of this inverse square root law in both cases of ordinary and fractional derivative formulations. Experimental affirmation of this law in optical bistable models within ordinary derivative formation was reported in [42].
In general, fractional order models provide deeper insight into the system dynamics with memory taken, into effect and further motivate for experimental observation.Finally, we refer to some very recent works [43,44] on biological models of COVID-19. In [43], the authors investigated various parameter estimation methods of COVID-19 incubation period using lognormal and Gamma distribution assumptions. The expressions for the maximum likelihood estimation, expectation maximisation algorithm and newly proposed algorithm [43] are termed as double or single (Riemann) integrals: these integral expressions can be converted to fractional integrals (i.e usual Riemann integral with memory or non-local, convolution kernel of fractional index, e.g. [23]), and so to have extra fractional order parameter. The other biological model of COVID-19 [44] is concerned with the stability and sensitivity analysis, and optimal control strategies of a suggested epidemic control of COVID-19. The adopted model of ODEs can be converted to FDEs and so to investigate the memory effect in this epidemic model. The formulation of the models in [43,44] within fractional calculus will certainly add details concerning memory/nonlocal effects.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.