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.

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 [15]. In general, researchers build models based on fractional differential by analogy to the approach that would be used in the case of classical differential [68]. 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 time t, 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 [1217] 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 function , its Riemann–Liouville fractional integral of order is defined aswith the provision that the right-hand side of the integral is point wise defined on and representing the usual gamma function .

Definition 2. Given a function , its Riemann–Liouville fractional derivative of order is defined bywhere .

Definition 3. Considering a function , its Caputo derivative of order is defined aswhere is the integer part of

Definition 4. Let be a continuous function and , then the Grünwald– Letnikov (GL) fractional derivative of is given bywhere

Definition 5. Let , then the new Caputo version of fractional derivative is defined aswhere is the normalization function with . If , then the new derivative called Caputo–Fabrizio fractional derivative can be defined as

Definition 6. Let . Then, the Atangana–Baleanu operator in the Caputo sense (ABC derivative) is defined aswhere is the normalization function with .

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 bywhere the fractional order of derivative and . Moreover, let us assume that the following functions , and are continuous on the closed interval . In order to compute the numerical solution to the problem defined by equation (9) over the interval , a discretization of into k subintervals of equal width is required. The set of points is also used in the approximation process. The GEM approximate solution is defined bywhere the node .

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 of derivative at the points , in the Grünwald–Letnikov sense has the following form [19]:where is the memory length; ; the time-space step of iteration is h and are referred to as binomial coefficients. For computational issue, the binomial coefficients are usually denoted by and computed as follows: .
Consider a nonlinear fractional differential equation, where the fractional derivative is taken in the Grünwald–Letnikov sense, with initial condition, defined byThe numerical solution to the problem stated by equation (12) is given by

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

The numerical solution to the problem defined by equation (14) is built based on the Adam–Bashforth method as (see [20, 21])

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

Let us consider the following fractional differential equation:

The numerical solution to the problem defined by equation (16) is built based on the Adam–Bashforth method as (see [20, 21])

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 equationswhere

The system given by equation (18) is converted into a fractional differential system by introducing the Caputo–Fabrizio fractional derivative as follows:

Alongside with the initial conditions,

Let us assume that the three variables involved in the process are such that , let be the Banach space of continuous functions defined on the interval endowed with the norm , where , , . Particularly, , where is the Banach space of continuous functions defined on .

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

It follows from equation (22) that

Let us denote for convenience,provided that the functions x, y, and z are bounded as assumed above.

Theorem 1. The kernel satisfies the Lipschitz condition and contraction, if the following inequality holds:

Proof. Let be two functions, thenHence, Lipschitz condition satisfied for and if additionally then it is also contraction.

Theorem 2. The kernel satisfies the Lipschitz condition and contraction, if the following inequality holds:

Theorem 3. The kernel satisfies the Lipschitz condition and contraction, if the following inequality holds:

Remark 1. Theorems 2 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:Let us consider the following system of recursive formula derived from equation (29):The difference between two consecutive terms of each of the equation from the system above can be written as follows:It follows from equations (30) and (31) thatMoreover, considering the first equation of the system defined by equation (31), we have the relationSince fulfills the Lipschitz condition, it follows from equation (33) thatSimilarly, it follows that

Theorem 4. The system of fractional differential equations for bone remodeling process equation (20) has an exact coupled solution, if there is such that

Proof. The following conditions hold on functions x, y, and z, , , and ; 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: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 thatIt follows from the assumption thatApplying the process recursively yieldsAt , we haveAt the limit case when n approaches infinity, it follows from equation (42) that .
Similarly, we can prove that and . 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 , thenApplying the norm on equation (43), it follows thatConsequently,

Theorem 5. The system of equation (20) has a unique solution if the following condition holds:

Proof. If the above condition holds, thenimplying thatSimilarly, we can prove that and .

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 .

The constants are ; ; .

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.

Figures 1(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 time t; 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.

Figures 2(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.

Figures 2(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 Section 5.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.

Figures 3(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.

Data Availability

No datasets were used or analyzed during the current study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Each of the authors, M.A, Y.Y.Y, and KA, contributed to each part of this work equally and read and approved the final version of the manuscript.

Acknowledgments

This article was funded by the Deanship of Scientific Research (DSR) under NASHER track with a grant number (206115), King Faisal University (KFU), Ahsa, Saudi Arabia. The authors, therefore, acknowledge the technical and financial support of DSR at KFU.