# A Fractional Epidemiological Model for Bone Remodeling Process

Authors: Muath Awadalla; Yves Yannick Yameni Noupoue; Kinda Abuasbeh
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This article focuses on modeling bone formation process using a fractional differential approach, named bones remodeling process. The first goal of the work is to investigate existence and uniqueness of the proposed fractional differential model. The next goal is to investigate how similar is the proposed approach to the method based on system classical differential equations. The dynamical system of equations used is built upon three main parameters. These are chemical substances, namely, calcitonin secretion, osteoclastic and osteoblastic, which are involved in the bone’s formation process. We implement some numerical simulations to graphically show the impact of an arbitrary fractional order of derivative. We finally obtained that modeling bone formation process using fractional differential equations yielded comparable results with those obtained through a system of classical differential equations. Flexibility in the choice of the fractional order of derivative is an advantage as it helps in selecting the best fractional order of derivative.


---

## Body

## 1. Introduction

Modeling natural phenomena through differential equation has long been used by scientists. In the early days, differential equations with integer order of derivative were commonly used. Then several fractional order derivatives have been implemented. In earlier days of fractional differential equations, researchers mainly focused on theoretical concepts, investigating existence and uniqueness of solution of built models. A sample of such theoretical works is found in [1–5]. In general, researchers build models based on fractional differential by analogy to the approach that would be used in the case of classical differential [6–8]. Hence, in many cases, authors will try to compare models’ performances from both approaches. Despite analogies that might exist between classical and fractional models upon a given problem, both models might not have common properties. For instance, the Carleman embedding technique [9] applicable in the studies of classical differential equations does not hold for the fractional differential equations for instance. Beside theoretical study of fractional differential equations, a lot have been done by researchers concerning possible application to real life phenomena. One can mention a nonexhaustive field of study with examples for which fractional differential has been successfully performed. In biology, modeling dynamics in human tissues was proposed in [10]. Alcohol can be dangerous when its concentration in high human blood exceeds a certain threshold. In [11], a fractional-based model of alcohol level in blood was proposed.Following the aims to prove the use of fractional differential equations, we decided to investigate how efficient it will be in modeling bone formation process. Before going into details, we will provide an overview of human bone formation process. Human bone formation and development are a complex process that starts from when a fetus is 3 months old and ends during teenage years, 13–18 years old. However, bone formation never ends in practice. Indeed, bones are living tissue made up of protein, calcium, and other minerals, as well as water. In this regard, bones’ tissue constantly renews itself, by breaking down older tissue and replacing it with new tissue. This process is called “bone remodeling process.” In literature, “bone formation” and “bone remodeling” refer often to the same process. There exit various cells or/and chemicals which are at the base of remodeling process. Three chemicals to which we will pay attention include the level of calcitonin above the basal level in the blood, the number of active osteoclasts at timet, and the number of active osteoblasts at time t.In this work, we build a system of fractional differential equations to modeling bone formation, which we call the bone remodeling process. Detailed information on the biological and chemical processes involved in the study as well as the chemical elements involved in bone formation are found in [12–17] and references therein. Therefore, we provided no or a very shallow biological and chemical description of the phenomenon. We rather focus on building a system of fractional differential equations of the model in which a classical differential equation counterpart is found in [12]. Moreover, we proved the existence and uniqueness of the built model. Lastly, we exhibited the solution of the built system using four numerical approaches of solving nonlinear differential equations, namely, the generalized Euler’s method (GEM), the Grünwald–Letnikov method or power series expansion (GL or PSE), the Caputo–Fabrizio method (CF), and the Atangana–Baleanu fractional derivative in Caputo sense (ABC). Modeling approach that consists of performing an experiment on the same data set with different techniques, including the proposed one, then to compare their performance using an error-metric function is not used in this work. Indeed, in the setting of this work, there is no comprehensive way to compare performances of the proposed method with those of existing methods.


## 2. Preliminaries

Fundamental definitions, terminologies, and notations commonly used in fractional calculus are provided in this section.Definition 1.
Given a functionh:0,+∞⟶ℝ, its Riemann–Liouville fractional integral of order q>0 is defined as(1)IRL0+qht=1Γq∫0tt−sq−1hsds,with the provision that the right-hand side of the integral is point wise defined on 0,+∞ and Γ representing the usual gamma function Γυ=∫0∞e−ttυ−1dt,∀υ>0.Definition 2.
Given a functionh:0,+∞⟶ℝ, its Riemann–Liouville fractional derivative of order q>0 is defined by(2)DRL0+qht=1Γn−qdndtn∫0tt−sn−q−1hsds,where n−1≤q<n,n∈ℕ.Definition 3.
Considering a functionh:0,+∞⟶ℝ, its Caputo derivative of order q>0 is defined as(3)DC0+qht=∫0tt−sn−q−1hnsΓn−qds,n−1<q<n∈ℝ,hnt,q∈ℕ,where n=q+1,q is the integer part of q.Definition 4.
Letg be a continuous function and n=t−a/h, then the Grünwald– Letnikov (GL) fractional derivative of g is given by(4)DGLa+qgt=limh⟶01hq∑j=0t−a/h−1jqjgt−jh,where(5)qj=q!j!q−j!=Γq+1Γj+1Γq−j+1,q0=1.Definition 5.
Letg∈H1a,b,b>aandq∈0,1, then the new Caputo version of fractional derivative is defined as(6)DCF0qgt=Mq1−q∫atg′sexp−q1−qt−sds,where Mq is the normalization function with M0=M1=1. If g∉H1a,b, then the new derivative called Caputo–Fabrizio fractional derivative can be defined as(7)DCF0qgt=Mq1−q∫0tg′sexp−q1−qt−sds.Definition 6.
Letg∈H1a,b,b>aandq∈0,1. Then, the Atangana–Baleanu operator in the Caputo sense (ABC derivative) is defined as(8)DABC0qgt=Bq1−q∫atg′sEq−q1−qt−sqds,where Bq is the normalization function with B0=B1=1.


## 3. Numerical Techniques of Solving Nonlinear Differential Equations

We discussed in this section some general numerical methods that are often used for computing numerical solution to nonlinear fractional differential equations. These methods are applied on the fractional model built in sequel.


### 3.1. The Generalized Euler’s Method (GEM)

The GEM is discussed in this subsection. This method was first introduced by Odibat and Momani [18]. The method is derived from the well-known Euler’s method of solving classical differential equations. Consider the fractional order nonlinear differential equation defined by(9)DC0+qt=gt,xt,x0=0,where the fractional order of derivative q∈0,1 and t>0. Moreover, let us assume that the following functions xt, DC0+qxt and DC0+2qxt are continuous on the closed interval 0,T. In order to compute the numerical solution to the problem defined by equation (9) over the interval 0,T, a discretization of 0,T into k subintervals tj,tj+1 of equal width h=T/k is required. The set of points tj,xtj is also used in the approximation process. The GEM approximate solution is defined by(10)xtj+1=xtj+hqΓq+1gtj,xtj,j=0,1,…,k−1,where the node tj=jh,j=1,2,…,k.


### 3.2. The Grünwald–Letnikov Method GL or PSE

Grünwald provided a numerical approach of solving nonlinear differential equation. This approach is discussed in detail in [19].Definition 7.
The explicit fractional numerical approximation formula ofqth derivative at the points kh, k=1,2,… in the Grünwald–Letnikov sense has the following form [19]:(11)DGLtkqgk−Lm/ht≈1hq∑j=0k−1jqjgtk−j,where Lm is the memory length; tk=kh; the time-space step of iteration is h and −1jqj are referred to as binomial coefficients. For computational issue, the binomial coefficients are usually denoted by cjq,j=0,1,… and computed as follows: c0q=1cjq=1−1+q/jcj−1q.
Consider a nonlinear fractional differential equation, where the fractional derivative is taken in the Grünwald–Letnikov sense, with initial condition, defined by(12)DGLa+qut=gut,t.
The numerical solution to the problem stated by equation (12) is given by(13)utk=gutk,tkhq−∑j=1kcjqutk−j.


### 3.3. The Caputo–Fabrizio Method (CF)

Consider a nonlinear fractional differential equation with initial condition, where the fractional derivative is taken in the Caputo–Fabrizio sense, defined by(14)DCF0qut=gt,ut,u0=u0.The numerical solution to the problem defined by equation (14) is built based on the Adam–Bashforth method as (see [20, 21])(15)un+1=un+1−qMq+3qh2Mqgtn,un+1−qMq+qh2Mqgtn−1,un−1.


### 3.4. Atangana–Baleanu Fractional Derivative in Caputo Sense (ABC)

Let us consider the following fractional differential equation:(16)DABC0qut=gt,ut,u0=u0.The numerical solution to the problem defined by equation (16) is built based on the Adam–Bashforth method as (see [20, 21])(17)un+1=un+1−qABCq+qABCqh2htn+1qq−2htn+1q+1q+1−qABCqΓhh2htn+1qq−2htn+1q+1q+1gtn,un+q−1ABCq−qABCqΓhhhtn+1qq−htn+1q+1q+1+tq+1ABCqΓhhgtn−1,un−1.


## 3.1. The Generalized Euler’s Method (GEM)

The GEM is discussed in this subsection. This method was first introduced by Odibat and Momani [18]. The method is derived from the well-known Euler’s method of solving classical differential equations. Consider the fractional order nonlinear differential equation defined by(9)DC0+qt=gt,xt,x0=0,where the fractional order of derivative q∈0,1 and t>0. Moreover, let us assume that the following functions xt, DC0+qxt and DC0+2qxt are continuous on the closed interval 0,T. In order to compute the numerical solution to the problem defined by equation (9) over the interval 0,T, a discretization of 0,T into k subintervals tj,tj+1 of equal width h=T/k is required. The set of points tj,xtj is also used in the approximation process. The GEM approximate solution is defined by(10)xtj+1=xtj+hqΓq+1gtj,xtj,j=0,1,…,k−1,where the node tj=jh,j=1,2,…,k.


## 3.2. The Grünwald–Letnikov Method GL or PSE

Grünwald provided a numerical approach of solving nonlinear differential equation. This approach is discussed in detail in [19].Definition 7.
The explicit fractional numerical approximation formula ofqth derivative at the points kh, k=1,2,… in the Grünwald–Letnikov sense has the following form [19]:(11)DGLtkqgk−Lm/ht≈1hq∑j=0k−1jqjgtk−j,where Lm is the memory length; tk=kh; the time-space step of iteration is h and −1jqj are referred to as binomial coefficients. For computational issue, the binomial coefficients are usually denoted by cjq,j=0,1,… and computed as follows: c0q=1cjq=1−1+q/jcj−1q.
Consider a nonlinear fractional differential equation, where the fractional derivative is taken in the Grünwald–Letnikov sense, with initial condition, defined by(12)DGLa+qut=gut,t.
The numerical solution to the problem stated by equation (12) is given by(13)utk=gutk,tkhq−∑j=1kcjqutk−j.


## 3.3. The Caputo–Fabrizio Method (CF)

Consider a nonlinear fractional differential equation with initial condition, where the fractional derivative is taken in the Caputo–Fabrizio sense, defined by(14)DCF0qut=gt,ut,u0=u0.The numerical solution to the problem defined by equation (14) is built based on the Adam–Bashforth method as (see [20, 21])(15)un+1=un+1−qMq+3qh2Mqgtn,un+1−qMq+qh2Mqgtn−1,un−1.


## 3.4. Atangana–Baleanu Fractional Derivative in Caputo Sense (ABC)

Let us consider the following fractional differential equation:(16)DABC0qut=gt,ut,u0=u0.The numerical solution to the problem defined by equation (16) is built based on the Adam–Bashforth method as (see [20, 21])(17)un+1=un+1−qABCq+qABCqh2htn+1qq−2htn+1q+1q+1−qABCqΓhh2htn+1qq−2htn+1q+1q+1gtn,un+q−1ABCq−qABCqΓhhhtn+1qq−htn+1q+1q+1+tq+1ABCqΓhhgtn−1,un−1.


## 4. Bone Remodeling Process

In this section, the bone remodeling process is discussed. Before the fractional model, we recall the classical differential equations model described in [12]. The system is made of three nonlinear differential equations(18)dxdt=c1+c2ym1+y−d1x,dydt=c3−c4xm2+x2yz−d2y,dzdt=c5+c6xm3+xz−d3z,where(19)xt:the level of CT above the basal level in the blood,yt:the number of active osteoclasts at timet,zt:the number of active osteoblasts at timet.The system given by equation (18) is converted into a fractional differential system by introducing the Caputo–Fabrizio fractional derivative as follows:(20)D0CFtqx=c1+c2ym1+y−d1x,D0CFtqy=c3−c4xm2+x2yz−d2y,D0CFtqz=c5+c6xm3+xz−d3z.Alongside with the initial conditions,(21)x0=x0,y0=y0,z0=z0.Let us assume that the three variables involved in the process are such thatT=X+Y+Z, let Ψ be the Banach space of continuous functions defined on the interval J endowed with the norm x,y,z=x+y+z, where x=supxt,t∈J, y=supyt,t∈J, z=supzt,t∈J. Particularly, Ψ=CJ×CJ, where CJ is the Banach space of continuous functions defined on J.


### 4.1. Existence and Uniqueness of the System of Bone Remodeling Process of Fractional Order

The existence and uniqueness of a solution for the system of fractional differential equations defined by equation (20) are investigated in this section through the fixed-point theory.It follows from the application of the fractional integral operator introduced by Nieto and Losada [22], to equation (20), that(22)xt−x0=I0CFtqc1+c2ym1+y−d1x,yt−y0=I0CFtqc3−c4xm2+x2yz−d2y,zt−z0=I0CFtqc5+c6xm3+xz−d3z.It follows from equation (22) that(23)xt−x0=21−q2−qMqc1+c2ytm1+yt−d1xt+2q2−qMq∫0tc1+c2ysm1+ys−d1xsds,yt−y0=21−q2−qMqc3−c4xtm2+x2tytzt−d2yt+2q2−qMq∫0tc3−c4xsm2+x2syszs−d2ysds,zt−z0=21−q2−qMqc5+c6xtm3+xtzt−d3zt+2q2−qMq∫0tc5+c6xsm3+xszs−d3zsds.Let us denote for convenience,(24)Q1t,x=c1+c2ytm1+yt−d1xt,Q2t,y=c3−c4xtm2+x2tytzt−d2yt,Q3t,z=c5+c6xtm3+xtzt−d3zt,μ1=d1,μ2=c3−c4xm2+x2z−d2,μ3=c5+c6xm3+x−d3,provided that the functions x, y, and z are bounded as assumed above.Theorem 1.
The kernelQ1 satisfies the Lipschitz condition and contraction, if the following inequality holds:(25)0≤μ1<1.Proof.
Letxandx1 be two functions, then(26)Q1t,x−Q1t,x1=c1+c2ytm1+yt−d1xt−c1+c2ytm1+yt−d1x1t≤μ1x−x1.
Hence, Lipschitz condition satisfied forQ1 and if additionally 0≤μ1<1, then it is also contraction.Theorem 2.
The kernelQ2 satisfies the Lipschitz condition and contraction, if the following inequality holds:(27)0≤μ2<1.Theorem 3.
The kernelQ3 satisfies the Lipschitz condition and contraction, if the following inequality holds:(28)0≤μ3<1.Remark 1.
Theorems2 and 3 are proved in a similar approach to what is used in the proof of Theorem 1.
Using the notations introduced above, the system given by equation (23) is written in a simple form as follows:(29)xt=x0+21−q2−qMqQ1t,x+2q2−qMq∫0tQ1s,xds,yt=y0+21−q2−qMqQ2t,y+2q2−qMq∫0tQ2s,yds,zt=z0+21−q2−qMqQ3t,z+2q2−qMq∫0tQ3s,zds.
Let us consider the following system of recursive formula derived from equation (29):(30)xnt=21−q2−qMqQ1t,xn−1+2q2−qMq∫0tQ1s,xn−1ds,ynt=21−q2−qMqQ2t,yn−1+2q2−qMq∫0tQ2s,yn−1ds,znt=21−q2−qMqQ3t,zn−1+2q2−qMq∫0tQ3s,zn−1ds.
The difference between two consecutive terms of each of the equation from the system above can be written as follows:(31)ϕxnt=xnt−xn−1t=21−q2−qMqQ1t,xn−1−Q1t,xn−2+2q2−qMq∫0tQ1s,xn−1−Q1s,xn−2ds,ϕynt=ynt−yn−1t=21−q2−qMqQ2t,yn−1−Q2t,yn−2+2q2−qMq∫0tQ2s,yn−1−Q2s,yn−2ds,ϕzn=znt−zn−1t=21−q2−qMqQ3t,zn−1−Q3t,zn−2+2q2−qMq∫0tQ3s,zn−1−Q3s,zn−2ds.
It follows from equations (30) and (31) that(32)xnt=∑i=1nϕxit,ynt=∑i=1nϕyit,znt=∑i=1nϕzit.
Moreover, considering the first equation of the system defined by equation (31), we have the relation(33)ϕxnt=xnt−xn−1t≤21−q2−qMqQ1t,xn−1−Q1t,xn−2+2q2−qMq∫0tQ1s,xn−1−Q1s,xn−2ds.
SinceQ1 fulfills the Lipschitz condition, it follows from equation (33) that(34)ϕxnt≤21−q2−qMqμ1ϕxn−1t+2q2−qMqμ1∫0tϕxn−1sds.
Similarly, it follows that(35)ϕynt≤21−q2−qMqμ2ϕyn−1t+2q2−qMqμ2∫0tϕyn−1sds,(36)ϕznt≤21−q2−qMqμ3ϕzn−1t+2q2−qMqμ3∫0tϕzn−1sds.Theorem 4.
The system of fractional differential equations for bone remodeling process equation (20) has an exact coupled solution, if there is t0 such that(37)21−q2−qMqμ1+2q2−qMqμ1t0<1.Proof.
The following conditions hold on functionsx, y, and z, x≤M1<∞, y≤M2<∞, and z≤M3<∞; moreover, we have proven that the kernels satisfied the Lipschitz condition, hence on consideration of equations (34)–(36) and by using the recursive method, we can derive the following nonrecursive formulae:(38)ϕxnt≤xn021−q2−qMqμ1++2q2−qMqμ1tn,ϕynt≤yn021−q2−qMqμ2++2q2−qMqμ2tn,ϕznt≤zn021−q2−qMqμ3++2q2−qMqμ3tn.
The system given by equation (39) stands as the solution of equation (20). Moreover, each equation of this system is continuous. Hence, the existence and continuity are proven. Now, let us prove that equation (39) is a solution of the fractional system equation (20). Let us assume that(39)xt−x0=xnt−Txnt,yt−y0=ynt−Tynt,zt−z0=znt−Tznt.
It follows from the assumption that(40)Txnt=21−q2−qMqQ1t,x−Q1t,xn−1+2q2−qMq∫0tQ1s,x−Q1s,xn−1ds≤21−q2−qMqQ1t,x−Q1t,xn−1+2q2−qMq∫0tQ1s,x−Q1s,xn−1ds≤21−q2−qMqμ1x−xn−1+2q2−qMqμ1x−xn−1t.
Applying the process recursively yields(41)Txnt≤21−q2−qMq+2q2−qMqtn+1μ1n+1a.
Att=t0, we have(42)Txnt≤21−q2−qMq+2q2−qMqt0n+1μ1n+1a.
At the limit case whenn approaches infinity, it follows from equation (42) that Txnt⟶0.
Similarly, we can prove thatTynt⟶0 and Tznt⟶0. This proves the existence of a solution of equation (20), which furthermore is equation (38).
Let us now prove the uniqueness of a system of solution of system equation (20). Let us assume there exists another set of solution to the system equation (20); let us denote this by x1t,y1t,andz1t, then(43)xt−x1t=21−q2−qMqQ1t,x−Q1t,x1+2q2−qMq∫0tQ1s,x−Q1s,x1ds.
Applying the norm on equation (43), it follows that(44)xt−x1t≤21−q2−qMqμ1xt−x1t+2q2−qMqμ1txt−x1t.
Consequently,(45)xt−x1t1−21−q2−qMqμ1−2q2−qMqμ1t≤0.Theorem 5.
The system of equation (20) has a unique solution if the following condition holds:(46)1−21−q2−qMqμ1−2q2−qMqμ1t>0.Proof.
If the above condition holds, then(47)xt−x1t1−21−q2−qMqμ1−2q2−qMqμ1t≤0,implying that(48)xt−x1t=0,thusxt=x1t.
Similarly, we can prove thatyt=y1t and zt=z1t.


## 4.1. Existence and Uniqueness of the System of Bone Remodeling Process of Fractional Order

The existence and uniqueness of a solution for the system of fractional differential equations defined by equation (20) are investigated in this section through the fixed-point theory.It follows from the application of the fractional integral operator introduced by Nieto and Losada [22], to equation (20), that(22)xt−x0=I0CFtqc1+c2ym1+y−d1x,yt−y0=I0CFtqc3−c4xm2+x2yz−d2y,zt−z0=I0CFtqc5+c6xm3+xz−d3z.It follows from equation (22) that(23)xt−x0=21−q2−qMqc1+c2ytm1+yt−d1xt+2q2−qMq∫0tc1+c2ysm1+ys−d1xsds,yt−y0=21−q2−qMqc3−c4xtm2+x2tytzt−d2yt+2q2−qMq∫0tc3−c4xsm2+x2syszs−d2ysds,zt−z0=21−q2−qMqc5+c6xtm3+xtzt−d3zt+2q2−qMq∫0tc5+c6xsm3+xszs−d3zsds.Let us denote for convenience,(24)Q1t,x=c1+c2ytm1+yt−d1xt,Q2t,y=c3−c4xtm2+x2tytzt−d2yt,Q3t,z=c5+c6xtm3+xtzt−d3zt,μ1=d1,μ2=c3−c4xm2+x2z−d2,μ3=c5+c6xm3+x−d3,provided that the functions x, y, and z are bounded as assumed above.Theorem 1.
The kernelQ1 satisfies the Lipschitz condition and contraction, if the following inequality holds:(25)0≤μ1<1.Proof.
Letxandx1 be two functions, then(26)Q1t,x−Q1t,x1=c1+c2ytm1+yt−d1xt−c1+c2ytm1+yt−d1x1t≤μ1x−x1.
Hence, Lipschitz condition satisfied forQ1 and if additionally 0≤μ1<1, then it is also contraction.Theorem 2.
The kernelQ2 satisfies the Lipschitz condition and contraction, if the following inequality holds:(27)0≤μ2<1.Theorem 3.
The kernelQ3 satisfies the Lipschitz condition and contraction, if the following inequality holds:(28)0≤μ3<1.Remark 1.
Theorems2 and 3 are proved in a similar approach to what is used in the proof of Theorem 1.
Using the notations introduced above, the system given by equation (23) is written in a simple form as follows:(29)xt=x0+21−q2−qMqQ1t,x+2q2−qMq∫0tQ1s,xds,yt=y0+21−q2−qMqQ2t,y+2q2−qMq∫0tQ2s,yds,zt=z0+21−q2−qMqQ3t,z+2q2−qMq∫0tQ3s,zds.
Let us consider the following system of recursive formula derived from equation (29):(30)xnt=21−q2−qMqQ1t,xn−1+2q2−qMq∫0tQ1s,xn−1ds,ynt=21−q2−qMqQ2t,yn−1+2q2−qMq∫0tQ2s,yn−1ds,znt=21−q2−qMqQ3t,zn−1+2q2−qMq∫0tQ3s,zn−1ds.
The difference between two consecutive terms of each of the equation from the system above can be written as follows:(31)ϕxnt=xnt−xn−1t=21−q2−qMqQ1t,xn−1−Q1t,xn−2+2q2−qMq∫0tQ1s,xn−1−Q1s,xn−2ds,ϕynt=ynt−yn−1t=21−q2−qMqQ2t,yn−1−Q2t,yn−2+2q2−qMq∫0tQ2s,yn−1−Q2s,yn−2ds,ϕzn=znt−zn−1t=21−q2−qMqQ3t,zn−1−Q3t,zn−2+2q2−qMq∫0tQ3s,zn−1−Q3s,zn−2ds.
It follows from equations (30) and (31) that(32)xnt=∑i=1nϕxit,ynt=∑i=1nϕyit,znt=∑i=1nϕzit.
Moreover, considering the first equation of the system defined by equation (31), we have the relation(33)ϕxnt=xnt−xn−1t≤21−q2−qMqQ1t,xn−1−Q1t,xn−2+2q2−qMq∫0tQ1s,xn−1−Q1s,xn−2ds.
SinceQ1 fulfills the Lipschitz condition, it follows from equation (33) that(34)ϕxnt≤21−q2−qMqμ1ϕxn−1t+2q2−qMqμ1∫0tϕxn−1sds.
Similarly, it follows that(35)ϕynt≤21−q2−qMqμ2ϕyn−1t+2q2−qMqμ2∫0tϕyn−1sds,(36)ϕznt≤21−q2−qMqμ3ϕzn−1t+2q2−qMqμ3∫0tϕzn−1sds.Theorem 4.
The system of fractional differential equations for bone remodeling process equation (20) has an exact coupled solution, if there is t0 such that(37)21−q2−qMqμ1+2q2−qMqμ1t0<1.Proof.
The following conditions hold on functionsx, y, and z, x≤M1<∞, y≤M2<∞, and z≤M3<∞; moreover, we have proven that the kernels satisfied the Lipschitz condition, hence on consideration of equations (34)–(36) and by using the recursive method, we can derive the following nonrecursive formulae:(38)ϕxnt≤xn021−q2−qMqμ1++2q2−qMqμ1tn,ϕynt≤yn021−q2−qMqμ2++2q2−qMqμ2tn,ϕznt≤zn021−q2−qMqμ3++2q2−qMqμ3tn.
The system given by equation (39) stands as the solution of equation (20). Moreover, each equation of this system is continuous. Hence, the existence and continuity are proven. Now, let us prove that equation (39) is a solution of the fractional system equation (20). Let us assume that(39)xt−x0=xnt−Txnt,yt−y0=ynt−Tynt,zt−z0=znt−Tznt.
It follows from the assumption that(40)Txnt=21−q2−qMqQ1t,x−Q1t,xn−1+2q2−qMq∫0tQ1s,x−Q1s,xn−1ds≤21−q2−qMqQ1t,x−Q1t,xn−1+2q2−qMq∫0tQ1s,x−Q1s,xn−1ds≤21−q2−qMqμ1x−xn−1+2q2−qMqμ1x−xn−1t.
Applying the process recursively yields(41)Txnt≤21−q2−qMq+2q2−qMqtn+1μ1n+1a.
Att=t0, we have(42)Txnt≤21−q2−qMq+2q2−qMqt0n+1μ1n+1a.
At the limit case whenn approaches infinity, it follows from equation (42) that Txnt⟶0.
Similarly, we can prove thatTynt⟶0 and Tznt⟶0. This proves the existence of a solution of equation (20), which furthermore is equation (38).
Let us now prove the uniqueness of a system of solution of system equation (20). Let us assume there exists another set of solution to the system equation (20); let us denote this by x1t,y1t,andz1t, then(43)xt−x1t=21−q2−qMqQ1t,x−Q1t,x1+2q2−qMq∫0tQ1s,x−Q1s,x1ds.
Applying the norm on equation (43), it follows that(44)xt−x1t≤21−q2−qMqμ1xt−x1t+2q2−qMqμ1txt−x1t.
Consequently,(45)xt−x1t1−21−q2−qMqμ1−2q2−qMqμ1t≤0.Theorem 5.
The system of equation (20) has a unique solution if the following condition holds:(46)1−21−q2−qMqμ1−2q2−qMqμ1t>0.Proof.
If the above condition holds, then(47)xt−x1t1−21−q2−qMqμ1−2q2−qMqμ1t≤0,implying that(48)xt−x1t=0,thusxt=x1t.
Similarly, we can prove thatyt=y1t and zt=z1t.


## 5. Experimental Studies and Numerical Simulations

In this section, we performed the numerical study of the solution of the systems of classical differential equation (18) of the bone formation process alongside the solution of the fractional differential equation (20) that we introduced for the same process. Computation and plotting are all done using mathematical software called MATLAB. The values of initials parameters and constants appearing in the equations are often determined experimentally by a biologist. For convenience, we used in our simulations the values retrieved from [12]. The initial conditions of the problem are defined as x0=1,y0=5,andz0=5.The constants area1=0.1;a2=0.5;a3=0.4;a4=0.7; a5=0.7;a6=0.085;k1=3;k2=5; k3=2;b1=0.1;b2=0.2;andb3=0.2.


### 5.1. Classical Model

In this section, solution to the problem defined by equation (18) is discussed. First order derivative equations are used in the system. Recall that the model has four variables, namely x, y, z, and t, where x, y, and z are the three chemicals involved in the process, respectively, the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. The parameter t represents the time involved in the bone’s formation process. 4-D representation would be appropriate tool to simultaneously visualize the behavior of x, y, and z with respect to t. However, only 1-D, 2-D, and 3-D representations are efficient for representation. In this regard, we choose a 2-D representation to show the variation of each of the substances x, y, and z, involved in the process as a function of time t. Figures 1(a)–1(c) represent the variation of each of the elements x, y, and z over time, whereas Figure 1(d) shows a joint and simultaneous variation of x and z.Figure 1
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through classical differential equations.
(a)(b)(c)(d)Figures1(a)–1(d) represent the dynamic of chemicals involved in bone formation process. Figures 1(a)–1(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. One can observe from Figures 1(a)–1(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and tends to 0 with time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 1(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.


### 5.2. Letnikov Fractional Model

Likewise the classical model case, we use a 2-D visualization to display the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given timet; and the number of active osteoblasts at a given time t (namely, x, y, and z, respectively). A numerical approach based of the fractional Letnikov model is used to produce variation graphs, which are shown in sequel. The initial parameters and constants are same with those used in the previous section for the classical approach, except the fractional order of derivative, which we selected as q = 0.95 in this study. Figures 2(a)–2(c) represent the variation of each of the elements x, y, and z over time, whereas Figure 2(d) shows a joint and simultaneous variation of x and z.Figure 2
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through fractional model of bone remodeling process by Letnikov.
(a)(b)(c)(d)Figures2(a)–2(d) represent the dynamic of chemicals involved in bone formation process. Figures 2(a)–2(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t obtained from Letnikov fractional differential model. One can observe from Figures 2(a)–2(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and converges to a small but nonzero value over time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 1(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.Figures2(a)–2(c) have the same sinusoidal shape as their counterpart from classical model. They also decrease over time. However, maximum (respectively minimum) value of the graph in classical model is larger (respectively smaller) than those obtained in the Letnikov model. Figures 1(d) and 2(d) have same shape; they both show same trend.


### 5.3. Caputo–Fabrizio Fractional Model

In this section, numerical solution based on Caputo–Fabrizio fractional model is discussed. The initial parameters and constants are like those used in Section5.3. Figures 3(a)–3(c) represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t (namely, x, y, and z), whereas Figure 3(d) is a joint and simultaneous variation of the variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at a given time t.Figure 3
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through fractional model of bone remodeling process by Caputo and Fabrizio.
(a)(b)(c)(d)Figures3(a)–3(c) represent the dynamic of chemicals involved in bone formation process. Figures 3(a)–3(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. It is shown from Figures 3(a)–3(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and tends to 0 with time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 3(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.


## 5.1. Classical Model

In this section, solution to the problem defined by equation (18) is discussed. First order derivative equations are used in the system. Recall that the model has four variables, namely x, y, z, and t, where x, y, and z are the three chemicals involved in the process, respectively, the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. The parameter t represents the time involved in the bone’s formation process. 4-D representation would be appropriate tool to simultaneously visualize the behavior of x, y, and z with respect to t. However, only 1-D, 2-D, and 3-D representations are efficient for representation. In this regard, we choose a 2-D representation to show the variation of each of the substances x, y, and z, involved in the process as a function of time t. Figures 1(a)–1(c) represent the variation of each of the elements x, y, and z over time, whereas Figure 1(d) shows a joint and simultaneous variation of x and z.Figure 1
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through classical differential equations.
(a)(b)(c)(d)Figures1(a)–1(d) represent the dynamic of chemicals involved in bone formation process. Figures 1(a)–1(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. One can observe from Figures 1(a)–1(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and tends to 0 with time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 1(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.


## 5.2. Letnikov Fractional Model

Likewise the classical model case, we use a 2-D visualization to display the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given timet; and the number of active osteoblasts at a given time t (namely, x, y, and z, respectively). A numerical approach based of the fractional Letnikov model is used to produce variation graphs, which are shown in sequel. The initial parameters and constants are same with those used in the previous section for the classical approach, except the fractional order of derivative, which we selected as q = 0.95 in this study. Figures 2(a)–2(c) represent the variation of each of the elements x, y, and z over time, whereas Figure 2(d) shows a joint and simultaneous variation of x and z.Figure 2
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through fractional model of bone remodeling process by Letnikov.
(a)(b)(c)(d)Figures2(a)–2(d) represent the dynamic of chemicals involved in bone formation process. Figures 2(a)–2(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t obtained from Letnikov fractional differential model. One can observe from Figures 2(a)–2(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and converges to a small but nonzero value over time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 1(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.Figures2(a)–2(c) have the same sinusoidal shape as their counterpart from classical model. They also decrease over time. However, maximum (respectively minimum) value of the graph in classical model is larger (respectively smaller) than those obtained in the Letnikov model. Figures 1(d) and 2(d) have same shape; they both show same trend.


## 5.3. Caputo–Fabrizio Fractional Model

In this section, numerical solution based on Caputo–Fabrizio fractional model is discussed. The initial parameters and constants are like those used in Section5.3. Figures 3(a)–3(c) represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t (namely, x, y, and z), whereas Figure 3(d) is a joint and simultaneous variation of the variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at a given time t.Figure 3
(a)–(d)) Dynamic of chemicals involved in bone formation process obtained through fractional model of bone remodeling process by Caputo and Fabrizio.
(a)(b)(c)(d)Figures3(a)–3(c) represent the dynamic of chemicals involved in bone formation process. Figures 3(a)–3(c), respectively, represent the variation over time of the level of CT above the basal level in the blood; the number of active osteoclasts at a given time t; and the number of active osteoblasts at a given time t. It is shown from Figures 3(a)–3(c) that chemical production decreases over time for all the three chemicals involved in the process. Moreover, the variation in chemical production shows a sinusoidal shape reflecting the high fluctuation in the process. The level of CT above basal level fluctuates on either side of the value 2.5 and tends to converge toward 2.5 over time. The number of active osteoclasts decreases very quickly and tends to 0 with time. Finally, the number of active osteoblasts fluctuates and tends to 1 with time. Figure 3(d) aims to investigate if there exists a correlation between any two of the chemicals produced in the process. Without loss of generality, correlation between variation over time of the level of CT above the basal level in the blood and the number of active osteoblasts at time t is investigated. The outcome of the investigation is that there is not a correlation in the way the produced chemical evolved over time.


## 6. Conclusion and Comments

The goal of this paper was to build a system of differential equations for bone formation process based on various fractional differential equations. Indeed, Letnikov and Caputo–Fabrizio fractional derivatives were used alongside classical differential equation method for comparison purposes. The work does not follow traditional quantitative modeling approach in which historical data exist, various models including a proposed model are used to fit the data; and their error rates are compared to pick the best model. We put all these models together rather to investigate if the proposed fractional differential equations-based models are effective. In the course of the work, we proved existence and uniqueness of the proposed system of differential equations; moreover, we studied the solution of the system using numerical approaches. From the numerical simulation of the classical model and those from the proposed fractional model, one can observe similarity in both approaches. This consolidates the results obtained in the theoretical part of the work. As a matter of fact, one can conclude that fractional differential equation is effective in modeling bone formation process. In further study, we will investigate advantages and disadvantages of using such method over classical approach.



---

*Source: 1614774-2021-11-10.xml*