Open access peer-reviewed chapter

Numerical Simulations of Some Real-Life Problems Governed by ODEs

Written By

N. H. Sweilam and T. A. Assiri

Submitted: October 26th, 2015 Reviewed: April 26th, 2016 Published: August 24th, 2016

DOI: 10.5772/63958

Chapter metrics overview

1,439 Chapter Downloads

View Full Metrics

Abstract

In this chapter, some real-life model problems that can be formulated as ordinary differential equations (ODEs) are introduced and numerically studied. These models are the variable-order fractional Hodgkin–Huxley model of neuronal excitation (VOFHHM) and other models with the variable-order fractional (VOF) time delay, such as the 4-year life cycle of a population of lemmings model, the enzyme kinetics with an inhibitor molecule model, and the Chen system model. A class of numerical methods is used to study the above-mentioned models such as non-standard finite difference (NSFD) and Adams-Bashforth-Moulton (ABM) methods. Numerical test examples are presented.

Keywords

  • Ordinary differential equations
  • Variable and fractional order real-life mathematical models
  • Adams-Bashforth-Moulton method
  • Non-standard finite difference method
  • Error analysis in numerical computations

1. Numerical simulations for VOFHHM of neuronal excitation

1.1. Introduction

Hodgkin and Huxley proposed the famous Hodgkin–Huxley equations which quantitatively describe the generation of action potential of squid giant axon [1].

CmdVmdt=IextIionic,E1

where Cm is the membrane capacitance, Vm is the intracellular potential, tis the time, Iionic is the net ionic current flowing across the membrane, and Iext is the externally applied current [2].

The ionic current described in Hodgkin and Huxley’s giant squid axon can be divided into three different cases: sodium current, potassium current, and leakage current, for example, voltage-gated persistent K+ current with four activation gates (resulting in the term n4, where n is the activation variable for K+), voltage-gated transient Na+ current with three activation gates and one inactivation gate (the term m3 ħ) [3].

The total ionic current is therefore

Iionic=GNa(VENa)INa+GK(VEK)IK+GL(VEL)IL,E2

where INa and IK are the currents through Na+ and K+ channels, respectively. The current IL is the leak current and denotes all the residue currents through a cell membrane other than Na+ and K+ currents. G denotes conductance value for each individual ionic component, and E is an equilibrium potential.

Hodgkin and Huxley defined the probability for a gate that can be found in its permissive state to depend on the membrane potential, therefore incorporating the voltage-dependent conductance [4, 5]. A probability pranging from 0 to 1 can be defined for each gate, but since a large population of channels (and therefore gates) has to be considered, this leads to the following:

1p(t)pαpβp(t)

where 1−p(t) is the fraction in non-permissive state, p(t) is the fraction in permissive state, and αpand βpare called rate functions.

Transitions between permissive and non-permissive states are assumed to obey the relation:

dpdt=αp(V)(1p)βp(V)p,p=m,n,,E3

where m and n are the activation variables for the sodium and the potassium, respectively, and represent the sodium de-inactivation.

In the literature, Alan Lloyd Hodgkin and Andrew Huxley authored a succession of five papers describing the nonlinear ODEs that model how action potentials can be initiated and propagated through an axon [15]. HHM can be considered one of the most successful mathematical models in quantity lively describing biological phenomena, and it is the method which can be used in deriving the model of a squid that is directly applicable to many kinds of neurons and other excitable cells [57].

Action potentials occur in excitable cells, including neurons, muscle cells, and endocrine cells. It is known that the brain is a complicated network of a tremendous number of neurons. A neuron has a very special shape, which is much different from usual sphere-shaped or disk-like cells [8]. A soma is the main body of neuron from which a long cable called an axon is extended. Neurons transmit and exchange electric signals called action potentials or spikes. Then, the electric signals or information is transmitted in the direction from a dendrite to an axon, see Figure 1 [8].

Figure 1.

Neurons cells.

Doi et al. [8] have briefly explained the framework of the Hodgkin–Huxley formalism to model the action potential generation of neuron sand of other excitable cells. In 2012, Siciliano [9] has studied different numerical methods to solve the HHM. Six different numerical methods are first introduced and compared using a simple and arbitrary ODE. In 2014, Sweilam and Nagy [6] have presented numerical method for solving fractional Hodgkin–Huxley model (FHHM), described as follows:

CDηV=IextG¯Nam3h(VENa)G¯Kn4h(VEK)GL(VEL),Dηp=αp(V)(1p)βp(V)p,p=m,n,,

where 0 < η ≤ 1.

In this work, the above-described model in [6] will be extended to VOFHHM, which is more general than fractional model.

The main aim of studying the VOFHHM is to show the behavior of the action potential when the derivative expressed as the VOF in order to explain the extent of this impact on the gating variables.

1.2. VOFHHM

It is well known that HHM is based on the parallel thought of a simple circuit with batteries, resistors, and capacitors [9]. Current can be carried through the circuit as ions passing through the membrane (resistors) or by charging the capacitors of the membrane. In this circuit, current flow across the membrane has two major components: one, IC(μA cm−2) associated with charging the membrane capacitance C(μF cm−2) and, the other, Iionic associated with the movement of specific types of ions across the membrane.

For the VOF model, we suggest that:

IC=CDη(t)V,E4

where Vis the membrane voltage ((mV), measured in volts), 0 < η(t) ≤ 1, and Dη(t) is defined by Grünwald–Letinkov approximation for the VOF derivative as follows:

Dη(t)y(t)=limh0hη(t)i=0[t]h(1)i(η(t)i)y(tih),E5

where [t] denotes the integer part of tand his the step size.

Substituting (4) into (1), we obtain:

CDη(t)V=IextIionic(t).E6

Transitions between permissive and non-permissive states are assumed to obey the relation:

Dη(t)p=αp(V)(1p)βp(V)p,p=m,n,,E7

The conductance for each current component can be written in the form [5]:

GNa=G¯Nam3h,GNa=G¯Kn4,GL=GL,E8

where g¯Naand g¯Kare the maximal conductance of sodium and potassium, respectively, and gLis a constant. Substituting (2) and (8) into (6), we obtain:

CDη(t)V=IextG¯Nam3h(VENa)G¯Kn4h(VEK)gL(VEL).E9

Eqs. (9) with (7) are a generalization of HHM [5], and it can be described as follows:

CDη(t)V=Iextg¯Nam3h(VENa)g¯Kn4h(VEK)gL(VEL),
Dη(t)p=αp(V)(1p)βp(V)p,p=m,n,.

The variable papproaches the steady-state value [5]:

p*(V)=αpαp+βp,p=m,n,E10

when the membrane potential Vremains constant.

1.3. Discretizations of the model using NSFD method

The aim of this part is to introduce numerical discretization for Eqs. (7) and (9). The NSFD method [1017] is used where the step size h in the FDM is replaced by a function ψ(h). Also, the Grünwald–Letinkov definition will be used with V(tk)=Vk, p(tk), p(tk) = pk, and then, we claim

i=0k=1ciη(t)pk+1i=αpVk(1pk)βpVkpk,p=m,n,,E11
Ci=0k+1ciη(t)Vk+1i=Img¯Namk12mk(VkENa)g¯Knk13nk(VkEK)gL(VkEL).E12

Doing some algebraic manipulation to Eqs. (11) and (12) yields the following equations:

pk+1=i=1k+1ciη(t)pk+1i+αpVk(1pk)βpVkpkc0η(t),p=m,n,E13
Vk+1=Iexti=1k+1ciη(t)Vk+1ig¯Namk12mk(VkENa)g¯Knk13nk(VkEK)gL(VkEL)Cc0η(t)E14

where c0η(t)=ψ(h)η(t),ψ(h)=eh1.

1.4. Numerical experiments

VOFHHM can be summarized neatly into four separate variable-order fractional ODEs with some supporting functions. These equations are described in Eqs. (7) and (9), where the rate functions are listed below [5]:

αn(V)=0.01(V+1)[exp(V+1010)]1,βn(V)=0.125exp(V80),
αm(V)=0.1(V+25)[exp(V+2510)]1,βm(V)=4exp(V18),
α(V)=0.07exp(V20),β(V)=[exp(V+3010)]1,

and

g¯Na=120mS/cm2,g¯K=36mS/cm2,gL=0.3mS/cm2,ENa=120mV,EK=12mV,EL=10.6mV,with the initial conditions V(0) = V0, and p(0) = p0 = p* (V0). The behavior of the neuron can be simulated for different initial values of V. Figure 2 shows an action potential simulated by the solver ode45, i.e., fourth-order Runge–Kutta method, and NSFD method at η(t) = 1, of the variable-order Hodgkin–Huxley equations for zero Iext, and V0 = −40 mV. Also, in this figure, all parts of the action potential are presented including the rapid upraise, downfall, and unexcitable phase, according to the number sequence as follows:

  1. Number 1 means that the action potential begins to fire, the sodium channels are now open, and an influx of sodium occurs.

  2. Number 2 means that the membrane becomes more positive at this point due to sodium entering the cell. Now, the potassium channels open and leave the cell.

  3. Number 3 means that the sodium channels have now become refractory and no more sodium is allowed to enter.

  4. Number 4 means that the sodium channel is still in its refractory stage and only the potassium ions are passing through the membrane. As the potassium ions continue to leave the cell, the membrane potential moves toward the resting value.

  5. Number 5 means that the potassium channels close and the sodium channels begin to leave the refractory phase and reset to its resting phase.

Figure 2.

An action potential created by using ode45 and NSFD method atη(t) = 1.

Finally, the potassium returns to its resting value where it can wait another action potential (for more details see [8, 9]). Moreover, in Figure 3, the gating variables are correctly constrained between 0 and 1 as expected using ode45 and NSFD method, respectively.

Figure 3.

Plot of the different gating variables using ode45 (left) and NSFD method (right) atη(t) = 1.

Figure 4.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) =0.99– (0.01/100)t.

Figures 48 show the approximated action potential and the gating variables using the NSFD method obtained for different values of (t), where η(t) = 0.99 − (0.01/100)t, η(t) = 0.87 − (0.01/100)t, η(t) =0.79 − (0.001)t, and η(t) = 0.67 − (0.001)t, respectively. These figures show an example of a numerically solved solution of the Hodgkin–Huxley equations. Figures 48 (left) are a waveform of the membrane potential. A pulsatile input applied at a time t = 5 ms ((ms) time measured in second) induces an action potential. Figures 47 (right) show the dynamics of all three gating variables m, n, and h. Sodium activation, m, changes much more rapidly than either h or n. Figure 9 (left) describes the relationship between the variables n and h, where n is the activation variable of potassium channel expressed in a dark color, while h is inactivation variable of sodium channel expressed in a light color. The action potential, n, be in stillness stage and begins to enter the cell until it reaches the highest activity, unlike the action potential, h, at the top of its activity. Step by step, the action potential, h, decreases its activities and continues to leave the cell and moves toward the resting value. At a certain moment in time, the intersection between them are occurred, where n increased while h decreased. Figure 9 (right) describes the relationship between the variables m and h, where m is the activation variable of potassium channel expressed in a dark color, while h is inactivation variable of sodium channel expressed in a light color. The dynamics variables m and h have the same dynamics similar to n and h, where m increased while h decreased. Figure 10 describes the dynamics of gating variables n and m in 3D. The variables n and m are the activation variables of K+ and Na+ ionic channels, respectively. They are expressed in a dark color and a light color, respectively. They are in stillness stage and begin to enter the cell until they reach the highest activity, where they are in the case of increasing.

Figure 5.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) =0.87– (0.01/100)t.

Figure 6.

The approximated action potential (left), and the gating variables (right) using the NSFD method at (t) =0.79– (0.001)t.

Figure 7.

The approximated action potential (left) and the gating variables (right) using the NSFD method at (t) =0.67– (0.001)t.

Figure 8.

An action potential created from the NSFD method using different values ofη(t).

Figure 9.

The waveforms of the gating variables n and h (left), and m and h (right) in 3D.

Figure 10.

The waveforms of the gating variables n and m in 3D.

Advertisement

2. Real-life models governed by nonlinear delay differential equations (DDEs)

2.1. Introduction

DDEs are differential equations in which the derivatives of some unknown functions at present time are dependent on the values of the functions at previous times. In real-world systems, delay is very often encountered in many practical systems, such as control systems [18], lasers, traffic models [19], metal cutting, epidemiology, neuroscience, population dynamics [20], and chemical kinetics [21]. Recent theoretical and computational advancements in DDEs reveal that DDEs are capable of generating rich and plausible dynamics with realistic parameter values. Naturally, occurrence of complex dynamics is often generated by well-formulated DDE models. This is simply due to the fact that a DDE operates on an infinite-dimensional space consisting of continuous functions that accommodate high-dimensional dynamics.

For example, the Lotka–Volterra predator prey model [22] with crowding effect does not produce sustainable oscillatory solutions that describe population cycles. However, the Nicholson’s blowflies model [23] can generate rich and complex dynamics. Delayed fractional differential equations (DFDEs) are correspondingly used to describe dynamical systems [23]. In recent years, DFDEs begin to arouse the attention of many researchers [7, 19, 20, 2527]. Simulating these equations is an important technique in the research, and accordingly, finding effective numerical methods for the DFDEs is a necessary process. Several methods based on Caputo or Riemann–Liouville definitions [28] have been proposed and analyzed. For instance, based on the predictor–corrector scheme, Diethelm et al. introduced the ABM algorithm [2931] and mean while some error analysis was presented to improve the numerical accuracy. In 2011, Bhalekar and Daftardar-Gejji [25] have been extended the ABM algorithm to solve DDEs of fractional order and presented numerical illustrations to demonstrate utility of the method. In 2012, Ma et al. [32] have presented the numerical solution of a VOF financial system which is calculated by using the ABM method.

This section presents class of numerical method for solving the variable-order fractional nonlinear delay differential equations (VOFDDEs). The main aim of this part is to study VOFDDEs numerically.

2.2. ABM method for the VOFDDEs

In this section, the ABM algorithm has been extended to study the VOFDDEs, where the derivative is defined in the Caputo VOF sense. Special attention is given to prove the error estimate of the proposed method. Numerical test examples are presented to demonstrate utility of the method. Chaotic behaviors are observed in variable-order one-dimensional delayed systems.

In the following, we apply the ABM predictor–corrector method to implement the numerical solution of VOFDDEs.

Let us consider the following VOF system:

Dtα(t)y(t)=f(t,y(t),y(tτ)),t[o,T],0<α(t)1.E15
y(t)=g(t),t[τ,0],E16

where fis in general a nonlinear function.

Also, consider a uniform grid {tn= nh: n= −k, −k+ 1,…, −1,0,1,…, N} where kand Nare integers

such that h= τ/k. Let

yh(tj)=g(tj),j=k,k+1,,1,0,E17

and note that

yh(tjτ)=yh(jhkh)=yh(tjk),j=0,1,2,,N.E18

Applying the integral Jtn+1α(tn+1), which is defined by the Riemann–Liouville variable-order integral as

Jxα(x)f(x)=1Γ(α(x))ax(xt)α(x)1f(t)dt,E19

on both sides of (15) and using (16), we claim to:

y(tn+1)=g(0)+1Γ(α(tn+1))0tn+1(tn+1ζ)α(tn+1)1f(ζ,y(ζ),y(ζτ))dζ.E20

Further, the integral in Eq. (20) is evaluated using product trapezoidal quadrature formula [2931, 33]. Then, we have the following corrector formula:

yh(tn+1)=g(0)+hα(tn+1)Γ(α(tn+1)+2)f(tn+1,yh(tn+1),yh(tn+1τ))+hα(tn+1)Γ(α(tn+1)+2)j=0naj,n+1f(tj,yh(tj),yh(tjτ),

or

yh(tn+1)=g(0)+hα(tn+1)Γ(α(tn+1)+2)f(tn+1,yh(tn+1),yh(tn+1k))+hα(tn+1)Γ(α(tn+1)+2)j=0naj,n+1f(tj,yh(tj),yh(tjk)E21

where

aj,n+1={nα(tn+1)+1(nα(tn+1))(n+1)α(tn+1),j=0,(nj+2)α(tn+1)+1+(nj)α(tn+1)+12(nj+1)α(tn+1)+1,1jn,1,j=n+1,}E22
yh(tjk)vn+1={δynm+1+(1δ)ynm+1,ifm>1,δyn+1p+(1δ)yn,ifm=1,}E23

0 ≤ δ< 1 and the unknown term yh(tn+1) appears on both sides of (21).

Due to nonlinearity off, Eq. (21) cannot be solved explicitly for yh(tn+1), so we replace the term yh(tn+1) on the right-hand side by an approximation yhp(tn+1),which is called predictor. The product rectangle rule is used in (20) to evaluate the predictor term

yhp(tn+1)=g(0)+1Γ(α(tn+1))j=0nbj,n+1f(tj,yh(tj),yh(tjτ),

or

yhp(tn+1)=g(0)+1Γ(α(tn+1))j=0nbj,n+1f(tj,yh(tj),yh(tjk)E24

where bj,n+1=hα(tn+1)α(tn+1)((nj+1)α(tn+1)+(nj)α(tn+1).

2.3. Error analysis of the algorithm

In this subsection, we aim to introduce the following lemma, which will be used in the proof of main theorem.

Lemma 2.1. [31] (a) Let zC1[0, T], then

|0tn+1(tn+1t)α(tn+1)1z(t)dtj=0nbj,n+1z(tj)|1α(tn+1)ztn+1α(tn+1)h.E25

(b) Let zC2[0, T], then

|0tn+1(tn+1t)α(tn+1)1z(t)dtj=0n+1aj,n+1z(tj)|Cα(tn+1)Trz''tn+1α(tn+1)h.E26

Proof. To prove statement (a), by construction of the product rectangle formula, we find that the quadrature error has the representation [31]

0tn+1(tn+1t)α(tn+1)1z(t)dtj=0nbj,n+1z(tj)=j=0njh(j+1)h(tn+1t)α(tn+1)1(z(t)z(tj))dt.

Apply the mean value theorem of differential calculus to the second factor of the integrand on the right-hand side of the above equation

|0tn+1(tn+1t)α(tn+1)-1z(t)dtj=0nbj,n+1z(tj)|
zj=0njh(j+1)h(tn+1t)α(tn+1)1(tjh)dt,
=zh1+α(tn+1)α(tn+1)j=0n(11+α(tn+1)[(n+1j)1+α(tn+1)
(nj)1+α(tn+1)](nj)α(tn+1),
=zh1+α(tn+1)α(tn+1)((n+1)1+α(tn+1)1+α(tn+1)j=0njα(tn+1)),
=zh1+α(tn+1)α(tn+1)(otn+1jα(tn+1)dtj=0njα(tn+1)).

Here, the term in parentheses is simply the remainder of the standard rectangle quadrature formula, applied to the function tα(tn+1)and taken over the interval [0, n + 1]. Since the integrand is monotonic, we may apply some standard results from quadrature theory to find that this term is bounded by the total variation of the integrand, thus

|0tn+1(tn+1t)α(tn+1)1z(t)dtj=0nbj,n+1z(tj)|1α(tn+1)z(n+1)α(tn+1)h1+α(tn+1).

Similarly, we can prove (b).

Now, let us consider that f(⋅) in (15) satisfies the following Lipschitz conditions with respect to its variables

|f(t,y1,u)f(t,y2,u)|L1|y1y2|,
|f(t,y,u1)f(t,y,u2)|L2|u1u2|,E27

where L1 and L2 are positive constants.

Theorem 1. Suppose the solution y(t) ∈ C2 [0, T], of the Eqs. (15) and (16), satisfies the following two conditions:

|0tn+1(tn+1t)α(tn+1)1Dtα(tn+1)y(t)dtj=0nbj,n+1Dtα(tn+1)y(tj)|Ctn+1γ1hδ1E28

|0tn+1(tn+1t)α(tn+1)1Dtα(tn+1)y(t)dtj=0naj,n+1Dtα(tn+1)y(tj)|Ctn+1γ2hδ2,

with some γ1, γ2 ≥ 0, and δ1, δ2 < 0, then for some suitable T> 0, we have

max0jN|y(tj)yh(tj)|=O(hq),

where q= min (δ1+ α(t), δ2), N= [T/h], and Cis a positive constant.

Proof. We will use the mathematical induction to prove the result. Suppose that the conclusion is true for j = 0, 1,…, n, then we have to prove that the inequality also holds for j = n + 1. To do this, we first consider the error of the predictor yn+1p. From (27), we have

|f(tj,y(tj),y(tjτ))f(tj,yj,υj)||f(tj,y(tj),y(tjτ))f(tj,yj,y(tjτ))|+|f(tj,yj,y(tjτ))f(tj,yj,υj)|L1hq+L2hq=(L1+L2)hq.

So

|y(tn+1)yn+1p|=1Γ(α(tn+1))|0tn+1(tn+1t)α(tn+1)1f(t,y(t),y(tτ))dtj=0nbj,n+1f(tj,yj,υj)|
1Γ(α(tn+1))|0tn+1(tn+1t)α(tn+1)1Dtα(tn+1)y(t)dtj=0nbj,n+1Dtα(tn+1)y(tj)|+1Γ(α(tn+1))j=0nbj,n+1|f(tj,y(tj),y(tjτ))f(tj,yj,y(tjτ))|+|f(tj,yj,y(tjτ))f(tj,yj,υj)|Ctn+1γ1Γ(α(tn+1))+1Γ(α(tn+1))j=0nbj,n+1(L1+L2)Chq,

and from

j=0nbj,n+1=j=0njh(j+1)h(tn+1t)α(tn+1)1dt=1α(tn+1)tn+1α(tn+1)1α(tn+1)Tα(tn+1),

we have

|y(tn+1yn+1p)|CTγ1Γ(α(tn+1))hδ1+CTα(tn+1)Γ(α(tn+1)+1)hq.

Since

j=0naj,n+1j=0nhα(tn+1)α(tn+1)(α(tn+1)+1)[(kj+2)α(tn+1)+1(kj+1)α(tn+1)+1(kj+1)α(tn+1)+1+(kj)α(tn+1)+1],
=0tn+1(tn+2t)α(tn+1)dt0tn+1(tn+1t)α(tn+1)dt,=1α(tn+1)0tn+1[(tn+1t)α(tn+1)](t)dt,=1α(tn+1)tn+1α(tn+1)Tα(tn+1),

we have

|y(tn+1)yn+1p|=1Γ(α(tn+1))|0tn+1(tn+1t)α(tn+1)1f(t,y(t),y(tτ))dtj=0naj,n+1f(tj,yj,υj)an+1,n+1f(tn+1,yn+1p,υn+1)|,
1Γ(α(tn+1))[|0tn+1(tn+1t)α(tn+1)1f(t,y(t),y(tτ))dtj=0naj,n+1f(tj,yj,y(tjτ))|+j=0naj,n+1|f(tj,yj,y(tjτ))f(tj,yj,υj)|+an+1,n+1|f(tn+1,y(tn+1),y(tn+1τ))f(tn+1,yn+1p,υn+1)|]1Γ(α(tn+1))[Ctn+1γ2hδ2+Chqj=0naj,n+1+an+1,n+1L(hq+CTγ1Γ(α(tn+1))hδ1+CTα(tn+1)Γ(α(tn+1)+1)hq)]
1Γ(α(tn+1))[CTγ2+Lhhα(tn+1)(α(tn+1)+1)+CLTγ1Γ(α(tn+1)+2)+CLTα(tn+1)hα(tn+1)α(tn+1)Γ(α(tn+1)+2)]hq=Chq.

2.4. Numerical examples

We present the purpose of this subsection to show that the proposed scheme designed provides good approximations for VOFDDEs. Throughout this subsection, we discuss four examples and their numerical solutions.

Example 2.1. Consider a VOFDDE:

Dtα(t)y(t)=2y(t2)1+y(t2)9.65y(t),
y(t)=0.5,t0E29

Figure 11.

The numerical behavior of system (29) and chaotic attractors atα= 1.

Figure 12.

The numerical behavior of system (29) and chaotic attractors atα(t) = 0.97.

Figure 13.

The numerical behavior of system (32) and chaotic attractors atα(t) =0.99– (0.01/100)t.

In Figures 11 and 1(a, b) show the solutions y(t) and y(t− 2) of the system (29), for α = 1 and h = 0.1, whereas Figure 1(c) shows phase portrait of the system, i.e., plot of y(t) versus y(t − 2) for the same value of α. In this figure, it observed that the system (29) shows aperiodic chaotic behavior. Moreover, Figure 12 shows the plot of the numerical solutions y(t) and the plot of y(t) versus y(t− 2), respectively, at α(t) = 0.97. In the following, we choose different cases for α(t). Figures 13 and 14 show the solutions y(t) and y(t − 2) of the given system for a(t) = 0.99 − (0.01/100)tand α(t) = 0.95 − (0.01/100)t, respectively. The chaotic portrait for these values of α(t) is shown. In Figure 15, we have decreased the value of α(t), where α(t) = 0.85 − (0.01/100)t, and observed that the system becomes periodic.

Figure 14.

The numerical behavior of system (29) and chaotic attractors atα(t) =0.95– (0.01/100)t.

Figure 15.

The numerical behavior of system (29) and chaotic attractors atα(t) =0.85– (0.01/100)t.

2.5. The four-year life cycle of a population of lemmings model

Lemmings are small rodents, usually found in or near the Arctic, in tundra biomes. It used to be that every three to 4 years, there were massive numbers of lemmings in the mountains and then the next year it was gone. Therefore, it has been an interesting thing to try to find out why. Are there factors that affect the lemming population such as climate, temperature, precipitation, and the like, that is, what we refer to as the lemming cycle. The modern scientific study of lemmings started with work carried out by the Norwegian Professor of Zoology Robert Collett, who at the end of the nineteenth century gathered a great deal of information about lemmings. To try to understand why lemmings fluctuate both regularly and extensively is indeed an important problem in ecology [34]. In [27], Tavernini has solved a model of the 4-year life cycle of a population of lemmings in an integer order. In [25], Bhalekar and Daftardar-Gejji have studied the model in a fractional order. In the following example [27], we study the extension of the lemming model:

Example 2.2. Consider the variable-order version of the four-year life cycle of a population of lemmings

Dtα(t)y(t)=3.5y(t)(1y(t0.74)19),y(0)19.00001,y(t)=19,t<0.E30

Figure 16.

The numerical behavior of system (33) and plots of the (t),y(t– 0.74) relation for α = 1.

Figure 17.

The numerical behavior of system (33), and plots of they(t),y(t– 0.74) relation forα(t) = 0.97.

Figure 18.

The numerical behavior of system (33) and the stretching phenomena forα(t) = 0.99 – (0.01/100)t.

Figure 16 shows the solutions y(t) and y(t− 0.74) of system (33) for α = 1 and h = 0.1 and shows phase portrait of y(t) versus y(t − 0.74) for the same value of α. Figure 17 shows the plot of the numerical solutions y(t) and the plot of y(t) versus y(t− 0.74), respectively, at (t) = 0.97. The numerical results for VOFDDEs at different values of α(t) are given in Figures 1822. These figures show the stretching phenomena between the numbers of lemmings in y(t) versus y(t − 2), for the values α(t) = 0.99− (0.01/100)t, α(t) = 0.98− (0.02/100)t; α(t) = 0.95− (0.01/100)t, α(t) = 0.87− (0.02/100)t, and α(t) = 0.85− (0.01/100)t, respectively. It is observed that the phase portrait gets stretched as the value of α(t) decreases, and this stretching is towards positive side of the axis.

Figure 19.

The stretching.

Figure 20.

The stretching phenomena forα(t) =0.98− (0.02/100)t. Phenomena forα(t) = 0.95 − (0.01/100)t.

Figure 21.

The stretching phenomena for α(t) = 0.87 − (0.02/100)t.

Figure 22.

The stretching phenomena for α(t) = 0.85 − (0.01/100)t.

2.6. The enzyme kinetics with an inhibitor molecule model

An enzyme inhibitor is a molecule that binds to an enzyme and decreases its activity or completely inhibits the enzyme catalytic activity. It is well known that all these inhibitors follow the same rule to interplay in enzyme reaction. Furthermore, there are many factors that affect enzyme’s activity, such as temperature and pH. Many drug molecules are enzyme inhibitors, so their discovery and improvement are an active area of research in biochemistry and pharmacology to protect enzyme from any change. Therefore, studying the enzyme kinetics and structure–function relationship is vital to understand the kinetics of enzyme inhibition that in turn is fundamental to the modern design of pharmaceuticals in industries [35]. The frequency conversion mechanism in enzymatic feedback systems has been investigated with computer simulations in 1984 by Okamoto and Hayashi [36]. In [25], Bhalekar and Daftardar-Gejji have been study the enzyme kinetics with an inhibitor molecule model in a fractional order. In the following, we will study the model in [25] in general form where the derivative is given in a VOF.

Example 2.3. Consider the variable-order version of four-dimensional enzyme kinetics with an inhibitor molecule

Dtα(t)y1(t)=10.5y1(t)1+0.0005y43(t4),
Dtα(t)y2(t)=y1(t)1+0.0005y43(t4)y2(t),
Dtα(t)y3(t)=y2(t)y3(t),
Dtα(t)y4(t)=y3(t)0.5y4(t),
y(t)=[60,10,10,20]T,t0.E31

Figure 23.

The numerical behavior of system (34) atα= 1 (left) andα= 95 (right).

Figure 23 shows the solutions y(t), (1≤ i≤ 4), with step size h= 1, for = 1, and α= 95, respectively. Whereas Figures 2426 show the numerical results in case of variable order at (t) = 0.99 − (0.01/100)t, α(t) = 0.95− (0.01/100)t, α(t) = 0.91− (0.01/100)t, and α(t) = 0.83− (0.01/100)t, respectively. For 0.90 < α(t) ≤ 1, the height of oscillations of the solutions increases as t increases as shown in Figures 24 and 25, while in Figure 26, the system settles down for sufficiently large t, for α(t) < 90.

Figure 24.

The numerical behavior of system (34) atα(t) = 0.99 – (0.01/100)t(left) andα(t) = 0.95 – (0.01/100)t(right).

Figure 25.

The numerical behavior of system (34) atα(t) = 0.91 – (0.01/100)t.

Figure 26.

The numerical behavior of system (34) atα(t) = 0.83 – (0.01/100)t.

2.7. The Chen system model

The atmosphere is a layer of gases surrounding the planet Earth that has, at each altitude above each point of the Earth’s surface, a density, a pressure, a temperature, etc., and all these vary over time. It is of course unthinkable to know all this infinite number of data and to understand something about it. Approximations must be made. Edward Norton Lorenz (1917–2008) is a pioneer of chaos theory. He introduced the strange attractor notion and coined the term butterfly effect. In 1963 [37], he developed a simplified mathematical model for atmospheric convection when he studied the atmospheric convection. His model of the atmosphere was reduced to only three numbers x, y, and z and the evolution of the atmosphere to a tiny equation, where each point x, y, and z in space symbolizes a state of the atmosphere and the evolution follows a vector field. Later, several dynamical systems exhibiting chaos have been presented in various branches of science [24]. For example, in 1999, Chen and Ueta found another simple three-dimensional autonomous system, which is not topologically equivalent to Lorenz’s system and which has a chaotic attractor [38]. In [39], Chen proved that “The Chen system is a special case of the Lorenz system.” In 2004, Li and Peng [40] have presented Chen system with a fractional order. In 2012, Bhalekar et al. [41] proposed the fractional-order Chen system with time delay. In the following example, we will extend Chen system to a VOF, which can be more general.

Example 2.4. Consider the generalization of the delay fractional-order version of the Chen system [3] which involves the variable order:

Dtα(t)x(t)=a(y(t)x(tτ)),
Dtα(t)y(t)=(ca)x(tτ)x(t)z(t)+cy(t),
Dtα(t)z(t)=x(t)y(t)bz(tτ),
x(t)=0.2,y(t)=0,z(t)=0.5fort[τ,0],E32

on the interval [0,30] and with step size h= 0.001.

At = 1, Figure 27 shows the numerical solution for τ= 0.005and shows chaotic yzphase portrait for this case. Figure 28 shows the numerical solution and chaotic yzphase portrait at α= 0.97and τ= 0.005. Moreover, for the variable order, we choose different cases for α(t):

At α(t) = 0.97– (0.01/100)t, Figure 29 shows the numerical solution y(t) and chaotic yzphase portrait of this example for τ= 0.005. When we increase the value of τ to 0.015, the results of the numerical solution and chaotic yzphase portrait become stable as shown in Figure 30. At α(t) = 0.94– (0.01/100)t, Figure 31 shows the numerical solution for τ= 0009, and a limit cycle phase portrait was observed. When we increase the value of τ to 0.011, the chaotic yzphase portrait becomes stable as shown in Figure 32. Figure 33 shows the numerical solution for τ= 0.005and chaotic yzphase portrait at α(t) = 89− (0.01/100)t. It is observed that even in one-dimensional delayed systems of variable order, chaotic behavior can be shown, and subjected to some critical order, the system changes its nature and becomes periodic. In some cases, it is observed that the phase portrait gets stretched as the order of the derivative is reduced. According to the numerical test examples, it can be concluded that the proposed numerical technique is a powerful technique to calculate approximate solution of VOFDDEs.

Figure 27.

The numerical behavior and chaotic attractors for =1, andτ=0.005.

Figure 28.

The numerical behavior and chaotic attractors forα = 0.97andτ = 0.005.

Figure 29.

The numerical behavior and chaotic attractors forα(t) = 0.97 – (0.01/100)tandτ = 0.005.

Figure 30.

The numerical behavior and chaotic attractors forα(t) =0.97– (0.01/100)tandτ=0.015.

Figure 31.

The numerical behavior and chaotic attractors forα(t) =0.94– (0.01/100)tandτ = 0.009.

Figure 32.

The numerical behavior and chaotic attractors forα(t) =0.94– (0.01/100)tandτ = 0.011.

Figure 33.

The numerical behavior and chaotic attractors forα(t) =0.89– (0.01/100)tandτ = 0.005.

References

  1. 1. A. L. Hodgkin, A. F. Huxley, B. Katz, Measurement of current-voltage relations in the membrane of the giant axon of Loligo, J. Physiol, 116, 4, 424–448, 1952.
  2. 2. A. L. Hodgkin, A. F. Huxley, Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo, J. Physiol, 116, 4, 449–472, 1952.
  3. 3. A. L. Hodgkin, A. F. Huxley, The components of membrane conductance in the giant axon of Loligo, J. Physiol, 116, 4, 473–496, 1952.
  4. 4. A. L. Hodgkin, A. F. Huxley, The dual effect of membrane potential on sodium conductance in the giant axon of Loligo, J. Physiol, 116, 4, 497–506, 1952.
  5. 5. A. L. Hodgkin, A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, J. Physiol, 117, 500–544, 1952.
  6. 6. N. H. Sweilam, A. M. Nagy, An efficient method for solving fractional Hodgkin-Huxley model, Phys. Lett. A, 378, 1980–1984, 2014.
  7. 7. N. H. Sweilam, T. A. Assiri, Numerical study for the general Hodgkin Huxley model of neuronal excitation, Appl. Math. Inf. Sci. 4, 8, 1–8, 2014.
  8. 8. S. Doi, J. Inoue, Z. Pan, K. Tsumoto, Computational electrophysiology, Springer, 2010.
  9. 9. Siciliano, Ryan. (2012). “The Hodgkin-Huxley Model – Its Extensions, Analysis and Numerics.” McGill University. Web. April 9th, 2013.
  10. 10. R. E. Mickens, Non-standard finite difference model of differential equations, World Scientific, Singapore, 1994.
  11. 11. R. E. Mickens, Application of non-standard finite difference schemes, World Scientific Publishing Co. Pte. Ltd., Singapore, 2000.
  12. 12. R. E. Mickens, Non-standard finite difference schemes for reactions-diffusion equations, Numerical Methods Partial Differential Equations Fractals, 15, 201–214, 1999.
  13. 13. R. E. Mickens, Numerical integration of population models satisfying conservation laws: NSFD methods, Biol. Dyn. 1, 4, 1751–1766, 2007.
  14. 14. R. E. Mickens, Calculation of denominator functions for non-standard finite difference schemes for differential equations satisfying a positivity condition, Numerical Methods for Partial Differential Equations, 23, 672–691, 2007.
  15. 15. R. E. Mickens, Finite difference models of ordinary differential equations: influence of denominator functions, J. Franklin Inst., 327, 143–149, 1990.
  16. 16. R. E. Mickens, A non-standard finite difference scheme for a fisher PDF having non-linear diffusion, Comput. Math. Appl., 45, 429–436, 2003.
  17. 17. R. E. Mickens, A non-standard finite difference scheme for the Lotka-Volterra system, Appl. Numer. Math., 45(2–3), 309–314, 2003.
  18. 18. E. Fridman, L. Fridman, E. Shustin, Steady modes in relay control systems with time delay and periodic disturbances, J. Dyn. Syst., Measur., Control, 122, 4, 732–737, 2000.
  19. 19. L. C. Davis, Modification of the optimal velocity traffic model to include delay due to driver reaction time, Physica A, 319, 557–567, 2002.
  20. 20. Y. Kuang, Delay differential equations with applications in population biology, Academic Press, Boston, San Diego, New York, 1993.
  21. 21. I. Epstein, Y. Luo, Differential delay equations in chemical kinetics, non-linear models: the cross-shaped phase diagram and the originator, J. Chem. Phys., 95, 244–254, 1991.
  22. 22. P. Liu, S. Elaydi, Discrete competitive and cooperative methods of Lotka-Volterra type, Comput. Appl. Analysis, 3, 53–73, 2001.
  23. 23. E. Braverman, D. Kinzebulatov, Nicholson’s blowflies equation with a distributed delay, Can. Appl. Math. Quart., 14, 2, 2006.
  24. 24. T. K. Alligood, T. D. Sauer, J. A. Yorke, Chaos: an introduction to dynamical systems, Springer, 2008.
  25. 25. S. Bhalekar, V. Daftardar-Gejji, A predictor-corrector scheme for solving non-linear delay differential equations of fractional order, J. Fract. Calculus Appl., 1, 5, 1–9, 2011.
  26. 26. K. S. Miller, B. Ross, An introduction to the fractional calculus and fractional differential equations, Wiley-Interscience, New York, NY, USA, 1993.
  27. 27. N. H. Sweilam, T. A. Assiri, A. M. Nagy, N. Y. Ali, Numerical simulations for variable order fractional non- linear delay differential equations, J. Fract. Calculus Appl., 6, 1, 71–82, 2015.
  28. 28. F. Liu, P. Zhuang, V. Anh, I. Turner, A fractional order implicit difference approximation for the space-time fractional diffusion equation, ANZIAM J., 47(EMAC2005), 48–68, 2006.
  29. 29. K. Diethelm, The analysis of fractional differential equations, Springer, Berlin, Germany, 2010.
  30. 30. K. Diethelm, N. J. Ford, A. D. Freed, A predictor-corrector approach for the numerical solution of fractional differential equations, Non-Linear Dynamics, 29, 3–22, 2002.
  31. 31. K. Diethelm, N. J. Ford, A. D. Freed, Detailed error analysis for a fractional Adams method, Numer. Algor., 36, 1, 31–52, 2004.
  32. 32. S. Ma, Y. Xu, W. Yue, Numerical solutions of a variable order fractional financial system, J. Appl. Math., 2012, 14, 2012 (Article ID 417942).
  33. 33. K. Diethelm, N. J. Ford, A. D. Freed, Yu. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Energy, 194, 743–773, 2005.
  34. 34. R. Boonstra, C. J. Krebs, N. C. Stenseth, Population cycles in small mammals: the problem of explaining the low phase, Ecology, 79, 1479–1488, 1998.
  35. 35. A. J. Sami, A. R. Shakoor, Cellulase activity inhibition and growth retardation of associated bacterial strains of Aulacophora foviecollis by two glycosylated flavonoids isolated from Magnifier indicia leaves, J. Med. Plants Res., 5, 2, 184–190, 2011.
  36. 36. M. Okamoto, K. Hayashi, Frequency conversion mechanism in enzymatic feedback systems, J. Theoret. Biol., 108, 4, 529–537, 1984.
  37. 37. N. E. Lorenz, Deterministic non-periodic flow, J. Atmos. Sci., 20, 1, 130–141, 1963.
  38. 38. G. Chen, T. Ueta, Yet another chaotic attractor, Int. J. Bifurc. Chaos, 9, 7, 1465–1466, 1999.
  39. 39. G. Chen, Comments on “Chen’s attractor exists if Lorenz repulse exists: The Chen system is a special case of the Lorenz system”, CHAOS, 23, 033108, 2013.
  40. 40. C. Li, G. Peng, Chaos in Chen’s system with a fractional order, Chaos, Solitons and Fractals, 22, 443–50, 2004.
  41. 41. S. Bhalekar, V. Daftardar-Gejji, P. Gade, Dynamics of fractional ordered Chen system with delay, PRAMANA-J. Phys., 79, 1, 61–69, 2012.

Written By

N. H. Sweilam and T. A. Assiri

Submitted: October 26th, 2015 Reviewed: April 26th, 2016 Published: August 24th, 2016