Numerical and graphical simulation of the non-linear fractional dynamical system of bone mineralization

1.   Introduction

Mathematical modeling is a powerful approach used in various scientific disciplines, including physics, engineering, economics, and biology ^[1]. It involves the construction of mathematical representations to describe, analyze, and predict real-world phenomena. In the context of biological systems, mathematical modeling plays a crucial role in obtaining an understanding of complex processes such as population dynamics, biochemical reactions, neural networks, and the spread of diseases.

In recent years, there has been a notable surge in the utilization of fractional calculus to model and elucidate biological processes. Fractional calculus has emerged as a sophisticated mathematical framework that enriches the traditional integer-order calculus by extending the concept of derivatives and integrals to non-integer orders (see, e.g., ^[2,3,4]).

The study of a dynamical system through the use of a fractional mathematical model provides a very convenient way of analyzing the phenomenon ^[5,6]. Modeling in fractional calculus helps us to completely understand the past course of the problem, evaluate the current situation, or predict the future course under different scenarios. Various mathematical models in fractional calculus have been studied by various researchers. For example, the dynamics of cytosolic calcium ions in astrocytes has been studied in ^[7,8], various models of epidemiology in ^[9,10], groundwater flow in karst aquifers in ^[11], modeling accuracy for complex systems like viscoelastic materials and anomalous diffusion processes in ^[12,13,14], and many others.

In 1967, Caputo ^[15] proposed a classic definition for a fractional derivative that overcomes the drawback of the Riemann-Liouville fractional derivative and allows the use of the classic initial conditions associated with integer-order differential equations. Caputo defined the fractional derivative in the following form:

Definition 1.1. Suppose that $\nu > 0, \, \, t > x, \, \, a, x \in \mathbb{R}$, and the fractional operator

is called the Caputo fractional derivative or Caputo fractional differential operator of order $\nu$, where $f(x) \in \mathrm{AC}^{n}[a, b]$.

Equivalently, in the convolved form, we have

Caputo and Fabrizio ^[16], in 2015, compiled available data regarding fractional derivatives, presenting a new branch of fractional calculus based on fractional operators with exponential, non-singular smooth kernels.

Let $H^1(a, b) = \{f|f \in L^2(a, b) \; \text{and}\; f'\in L^2(a, b)\}$, where $L^2(a, b)$ is the space of square-integrable functions on the interval $(a, b)$.

Definition 1.2. Let $\nu \in(0, 1)$, $f(x)\in H^{1}(a, b)$, and the Caputo-Fabrizio (C–F) fractional derivative be defined as

where $Z(\nu)$ is the normalization function such that $Z(0) = Z(1) = 1$.

The C–F operator, in terms of convolution, can be expressed as follows:

Losada and Nieto ^[17] derived the appropriate C–F fractional integral. The formal definition is as follows:

Definition 1.3. The C–F fractional integral for the function $f$ of order $\nu \in [0, 1)$ is given by

Fractional-order models differ from classical integer-order models in that they can address non-locality, anomalous diffusion, nonlinear phenomena, memory effects, lag, etc, in various science and engineering phenomena, taking into account the history and locally dispersed effects (see, e.g., ^[2,3]). In this study, we have attempted to modify the standard bone mineralization model by employing the C–F fractional differential operator.

The structure of bone material can be presented by the interaction of two dynamic processes known as mineralization and remodeling. Remodeling is the process whereby old, initially rigid bone is replaced with new, softer material. On the other side, mineralization refers to the occurrence of inorganic precipitation over an organic background. Since the bone material stiffens as a result of the mineralization process, this process plays an essential role in the stability of the bone ^[18].

Bone mineralization is a complex process that involves the deposition of minerals, primarily calcium and phosphate, onto the extracellular matrix of bone. The presence of calcium phosphate minerals in bone was initially noted by Scheele and Dobbin in 1771 ^[19]. The organic bone matrix crystallizes hydroxyapatite $[Ca_{10}(PO_4)_6(OH)_2]$ as a result of the precipitation of calcium and phosphate during the mineralization process ^[20]. The collagen matrix stiffens as a result of the mineralization process. The degree of mineralization is an essential determinant of bone quality since it modifies the matrix material's elastic modulus; the higher the degree of mineralization, the stiffer the material ^[21]. In addition to the organic matrix's material characteristics, the amount of mineral is an essential component in the stiffness of the bone material ^[18]. According to current theories, mineral clusters occur in the surrounding fluid and are thought to attach to specific regions of the collagen matrix as the first steps in the mineralization of newly deposited collagen. At the level of mineralized collagen fibrils, the elastic characteristics of bone are discussed in ^[22]. A novel modeling technique that enables the computation of physical characteristics of the human cortical bone is discussed in ^[23]. The mineralization profiles within a particular bone using various techniques are discussed in ^[24], and the development of dentin and dentinal tubules and their relation to mineralization and diffusion processes are covered in ^[25].

The dynamics of mineral deposition and resorption during bone mineralization can be viewed from a new perspective by using fractional calculus ^[26]. Motivated by the results obtained by researchers using fractional operators for various biological models, in this paper, we have fractionalized Komarova's bone mineralization model by using the C-F derivative and investigated the dynamics of the new model qualitatively.

This article is organized into various sections, as follows: Section 2 introduces the mathematical model of the dynamics of bone mineralization; in Section 3, the stability, existence, and uniqueness of the solution are discussed; Section 4 describes the numerical and graphical simulations and the sensitivity analysis with respect to the various parameters of the model; finally, Section 5 is the conclusion section.

2.   Mathematical model

Komarova et al. ^[27], in 2015, proposed the following mathematical bone mineralization model, which successfully describes the nonlinear mineralization dynamics. Komarova et al. described the dynamics of bone mineralization by applying integer order and explained the behavior of each variable. We consider the C–F derivative that has non-singular kernels and hence describes the distributed time delay and lag phenomenon ^[28]. The fractionalized mathematical model describes a process in a more realistic way than the integer order. Mineralization processes start with a lag and have a finite speed. The change in the number of raw collagen does not lead to instantaneous mineralization. The developed model is highly efficient in terms of predicting mineralization disorders. The model contains $5$ important variables, including raw collagen, mature collagen, inhibitors, nucleators, and mineralization.

The following system of equations describes the dynamics of bone mineralization:

The variables and parameters of the model are defined as follows:

$\phi$ is the concentration of the raw collagen,

$\psi$ is the concentration of the mature collagen,

$\theta$ is the number of inhibitors that prevent raw collagen from converting into mature collagen,

$\pi$ is the number of nucleators that act on the mature collagen so that mineralization occurs,

$\omega$ is the mineral,

$\kappa$ is the rate at which collagen cross-linking takes place and is inversely related to time lag during mineralization,

$\beta$ is the rate at which the inhibitors are released into the extracellular compartment near the cells and diffuse through immature collagen,

$\epsilon$ is the rate of removal of inhibitors and is stimulated in the presence of a mature collagen matrix,

$\lambda$ is the number of nucleators per mature collagen molecule,

$\upsilon$ is the rate of nucleator removal,

$\eta$ is the rate at which mineralization takes place.

The following presumptions are used in the mathematical modeling of the mineralization process:

● Matrix mineralization is modeled with assumed homogeneity.

● Inhibitors prevent raw collagen from converting into mature collagen.

● A nucleator is required to initiate mineralization, which happens quickly.

● Nucleators produced during collagen maturation are eliminated from the system at a rate that is proportional to the rate of mineralization.

Agarwal et al. ^[29] studied the integer-order model of bone mineralization. We hereby study the fractionalized model involving the C–F derivative, which has a non-singular kernel.

Upon replacing the integer-order derivative in (2.1) with C–F fractional derivative and balancing the dimensions ^[30], we obtain

To make the system dimensionless, we apply the following substitutions:

the system (2.2) is reduced as follows:

In the next sections, we analyze the model (2.4) qualitatively and numerically, including the boundedness, existence, and uniqueness of the solution and stability analysis of the coupled system of equations defined above. After that the numerical simulations and sensitivity analysis are presented.

3.   Qualitative analysis

This section examines the boundedness, existence, and uniqueness of the system of the differential equations (2.4)–(2.8). The contraction mapping theorem is used to establish the solution's existence. Subsequently, the stability analysis is carried out to determine the future behavior of the system ^[31].

3.1.   Boundedness of the solution

We must demonstrate that the model is biologically feasible and the parameter values are constrained before moving to the analytical approach.

Taking the Laplace transform of both sides of the equation (2.4), and using the formula

we obtain

Let

,

Hence, $\Phi(t)$ is bounded and positive.

Put the value of $\Phi$ from (3.4) in (2.5); then, we get

Taking the C–F integral, we get

By incorporating the values of $p, q$ from (3.3), we get

This implies that $\Psi(t)$ is bounded.

Given that $\hat{\epsilon}\Theta\Psi > 0$, from (2.6), we observe that ${}_0^{CF}D_t^{\nu}\Theta\leq\hat{\beta} \Phi$.

Applying the value of $\Phi$ from (3.4), we get

Applying the C–F integral, we get

Hence, $\Theta(t)$ is bounded.

Given that ${}_0^{CF}D_t^{\nu}\Omega = \hat{\eta} \bigg(\dfrac{\tilde{\rho}}{\tilde{\rho}+\Theta^\sigma}\bigg)\Pi > 0$, from (2.7), we observe that ${}_0^{CF}D_t^{\nu}\Pi\leq \hat{\kappa}\lambda\Phi$.

Applying the value of $\Phi$ from (3.4), we get

Taking the C–F integral, we get

Hence, $\Pi(t)$ is bounded.

Taking the C–F integral of (2.8), we get

since $\bigg(\dfrac{\tilde{\rho}}{\tilde{\rho}+\Theta^\sigma}\bigg) < 1$.

Applying the value of $\Pi(t)$ from (3.12) and integrating, we get

It can be observed that $\Omega(t)$ is bounded for $t$ lying in a bounded interval.

3.2.   Existence and uniqueness of the solution of the Bone mineralization system

In this subsection, we prove the existence and uniqueness of the solution of the nonlinear system of the differential equations (2.4)–(2.8).

If $C(J)$ is the collection of continuous real-valued functions defined on the interval $J \subset \mathrm{R}$. Then, $V = C(J)\times C(J) \times C(J)\times C(J) \times C(J)$ will be a Banach space with the norm for $(\Phi, \Psi, \Theta, \Pi, \Omega) \in V$, defined as $\|(\Phi, \Psi, \Theta, \Pi, \Omega)\| = \|\Phi\|+\|\Psi\|+\|\Theta\|+\|\Pi\|+\|\Omega\|$, where $\Phi, \Psi, \Theta, \Pi, \Omega \in C(J)$ and the norm is the supremum norm, which is defined as $\|\cdot\| = \sup_{t \in J}|\cdot|$.

Here, using the contraction mapping theorem, we prove the existence of the solution of the system of differential equations given by (2.4)–(2.8).

Applying the C–F integral operator to (2.4)–(2.8) respectively, we obtain

Upon using the definition of ${}^{CF}_0I_t^\nu$, we have

Let us write the following as kernels:

Theorem 3.1. The Lipschitz condition and contraction would be satisfied by $\mathfrak L_{1}, \mathfrak L_{2}, \mathfrak L_{3}, \mathfrak L_{4}, \mathfrak L_{5}$ if the inequalities $0\leq \hat{\kappa} < 1, \, \, 0\leq \dfrac{2\hat{\epsilon}\Phi_0}{Z(\nu)} < 1$, and $\, \, 0\leq 2\hat{\upsilon}\lambda \hat{\eta}\xi < 1$ hold.

Proof. Let us start with $\mathfrak L_{1}$, and with $\Phi$ and $\Phi^{(1)}$ as two functions satisfying (2.4) of the bone mineralization model. Then,

Clearly, $\hat{\kappa}$ is a fixed parameter and $\Phi(u)$ is a bounded function. Hence, it satisfies the Lipschitz condition for $\mathfrak L_{1}$, and it is contraction if $0 \leq\hat{\kappa} < 1$.

The Lipschitz condition is satisfied in the remaining four cases as well, i.e.,

□

Theorem 3.2. The system of equations to model fractional-order bone mineralization, as given by (2.4)–(2.8), has an exact coupled solution under the condition that there exists a $t_0$ such that $\dfrac{1-\nu}{Z(\nu)} \gamma_i+\dfrac{\nu}{Z(\nu)} \gamma_it_0 < 1$ holds for

and the solution is unique.

Proof. Existence:

Upon using the above kernels given in (3.21) in the equations (3.16)–(3.20) respectively, we get

The corresponding recursive formulas are respectively given by

The initial conditions are $\Phi^{(0)} = \Phi(0), \, \, \Psi^{(0)} = \Psi(0), \, \, \Theta^{(0)} = \Theta(0), \, \, \Pi^{(0)} = \Pi(0)$, and $\, \, \Omega^{(0)} = \Omega(0)$.

The following expression respectively represents the differences between the above terms in (3.26) and their succeeding terms:

Hence,

Further, we have

Upon applying the triangle inequality, we get

The kernels satisfy the Lipschitz condition given by

and, hence,

Similarly,

□

Since the functions $\Phi, \Psi, \Theta, \Pi, \Omega$ are bounded and the kernels satisfy the Lipschitz condition, from the equations (3.36)–(3.40) and the recursive method, the resulting relations can be derived as follows:

Therefore, the functions given by (3.27)–(3.31) will exist and be continuous.

Now, to show that the functions making up (3.32) constitute a solution of the system of the equations given by (2.4)–(2.8), consider the following relations:

Now,

Upon repeated application of the recursive relations, we get

At some $t_0$,

Taking the limit $n \to \infty$ in the equation (3.45), we have

Similarly,

This proves that the solution exists for the model.

Uniqueness:

For the uniqueness, let us consider $\{\Phi', \Psi', \Theta', \Pi', \Omega'\}$ as another solution of the system of equations given by (2.4)–(2.8):

It gives

Since $\bigg(1-\dfrac{1-\nu}{Z(\nu)} \gamma_1-\dfrac{\nu}{Z(\nu)} \gamma_1t\bigg) > 0,$

Similarly, we can show that $\Psi(t) = \Psi'(t), \; \Theta(t) = \Theta'(t), \; \Pi(t) = \Pi'(t)$, and $\; \Omega(t) = \Omega'(t).$

Hence, the uniqueness of the solution is proved.

3.3.   Stability analysis

A requirement of studying dynamic systems and understanding their behavior is stability analysis. The stability of the equilibrium points of the fractional-order differential system has been studied by using the Jacobian and corresponding eigenvalues at the points of equilibrium.

We discuss here the stability of Komarova's coupled system for the process of bone mineralization.

Theorem 3.3. The system given by (2.4)–(2.8) is marginally stable at the critical points.

Proof. In the model given by (2.4)–(2.8), we set all of the derivatives equal to zero and define them as follows:

Critical points can be obtained by setting $f_i$ = 0 for $i = 1, 2, 3, 4, 5$.

Since $\Phi+\Psi = K$, we have that $\Psi = K$.

Now,

Also,

${}_0^{CF}D_t^{\nu}\Omega$ may or may not be zero. Hence, the system has infinitely many critical points $(\Phi, \Psi, \Theta, \Pi, \Omega) = (0, K, 0, 0, \Omega)$.

We shall now determine whether the system in question is stable or unstable at its critical points. So, for this, we will find the Jacobian matrix.

For the particular system of bone mineralization, the Jacobian matrix will take the following form:

Making the substitutions for $f_i$, $i = 1, 2, 3, 4, 5$ in the matrix we get the following Jacobian matrix:

The eigenvalues corresponding to the above matrix are $(0, 0, 0, -\hat{\kappa}_1, -\hat{\kappa}_1\hat{\epsilon}_1)$.

Observing the eigenvalues, we can conclude that the system is marginally stable at the critical points $(0, K, 0, 0, \Omega)$, which occur after a short span of the start of the mineralization process. Thus, the mineralization will not suddenly explode and the system will always have a bounded solution but no steady-state output. □

4.   Numerical and graphical simulation

We solve the initial-value fractional-order problem by using the numerical technique. Here, we use the fractional Adams-Bashforth (AB) method ^[32] to find the approximation solution.

Let us consider the problem with the exponential law kernel:

For $Z(\nu) = \dfrac{2}{2-\nu}$ in the equation (1.5), the Lagrangian approximation of $w$, i.e., the solution of the equation (4.1), is given by

Following the process, the Lagrangian approximations for $\Phi, \Psi, \Theta, \Pi, \Omega$ in (2.4)–(2.8) are respectively given by

The following values were taken from ^[27] for the parameters before the non-dimensionalization:

$\phi = 9.4\times 10^5 \; \text{molecules}/\mu m^3, \; \; \psi = 9.4\times 10^5 \; \text{molecules}/\mu m^3, \; \; \theta \sim 10^6 \; \text{molecules}/\mu m^3$,

$\pi = 1-10 \, \, \text{per} \, \, 1\, \, \text{assembled collagen}, \; \; \omega = 0.8 \times10^9 \; \text{molecules}/\mu m^3$.

$\kappa = 0.1 \, \, \text{per day}, \; \; \lambda = 1, \quad \epsilon = 2 \times 10^{-7}\, \, \text{per day}, \; \; \eta = 1000\, \, \text{per day}, \; \; \beta = 0.1 \, \, \text{ per day},$

$\upsilon = 1.7\times 10^{-8}\, \, \text{per mol}, \; \; \sigma = 10, \; \; \rho = 10^{57}$.

For normalization, the following values have been considered in (2.3):

$\hat{\phi} = 10^6, \; \; \hat{\psi} = 10^9$.

Figure 1 shows how the concentrations of raw collagen, mature collagen, inhibitors, nucleators, and minerals change with respect to time. Here, the dotted line denotes $\nu = 1$ (integer order) and the solid line denotes $\nu = 0.9$ (fractional order). Initially, the raw collagen $(\Phi)$ has full collagen in the system; with time, it decreases gradually and is converted into mature collagen $(\Psi)$. For the orders $\nu = 0.9$ and $\nu = 1$, there is a little difference in the process. At the beginning, the inhibitors $(\Theta)$ are present in the raw collagen and then decrease rapidly within $20$ days. For the orders $\nu = 0.9$ and $\nu = 1$, there is a difference in the peak values of the inhibitors. For $\nu = 1$, the number of inhibitors is larger within $20$ days than for $\nu = 0.9$. The nucleators $(\Pi)$ reach the maximum within $20$ days and then decrease rapidly. For $\nu = 1$, it reaches the maximum peak faster than the case of $\nu = 0.9$. For the fractional order, the peak is small and it begins decreasing very slowly. The mineralization $(\Omega)$ starts with some lag time. For fractional order $\nu = 0.9$, the lag time is much shorter than for the integer order $\nu = 1$, which means that it takes less time for the mineralization. For the integer order $\nu = 1$ the degree of mineralization is a little higher than for the fractional order $\nu = 0.9$.

The amount of time needed to initiate the mineralization process is known as the Mineralization lag time, and the Mineralization degree is defined as the greatest amount of mineralization that may occur during the process.

The concentration of raw collagen with respect to time for the different values of the order $\nu$ of the fractional derivative is displayed in Figure 2a. Because the value of $\nu$ decreases during a specific period, raw collagen requires fewer days to be converted into mature collagen. The amount of mature collagen can be observed to change over time, as seen in Figure 2b. During a particular interval of time, the concentration of the mature collagen increases as $\nu$ decreases. The average life cycle for the conversion of raw collagen to mature collagen is less than 100 days; the curves in 2 merge after the completion of the cycle.

Figure 3a shows how inhibitors degrade with time for various values of $\nu$. The peak shifts and the number of inhibitors increases as the value of $\nu$ increases. For different values of $\nu$, the effect of nucleators can be observed in Figure 3b. The number of nucleators for mineralization increases as $\nu$ increases and the peak shifts in a small range. Figure 4 demonstrates the impact of mineralization for different values of $\nu$. As the value of $\nu$ increases the mineralization lag time and the degree of mineralization both increase. Initially, the mineralization lag time increases gradually, but when $\nu = 1$, it increases at a comparatively fast pace.

4.1.   Sensitivity analysis

Sensitivity analysis of a model with respect to its parameters determines how different values of an independent variable affect a particular dependent variable under a given set of assumptions. We present here the concentration profiles of the variables and the sensitivity analysis with respect to the parameters $\kappa, \beta, \epsilon, \lambda, \upsilon, \eta$ to illustrate their impact on the variables and, ultimately on the mineralization.

The effect of the parameter $\kappa$ on the raw collagen is shown in Figure 5a. It shows that, as the value of $\kappa$ increases, raw collagen begins to require fewer days to convert into mature collagen. The effect of the parameter $\kappa$ on the mature collagen is shown in Figure 5b. It demonstrates that, when the value of $\kappa$ increases, the time it takes for mature collagen to form decreases.

The effect of the parameter $\beta$ on the inhibitors is shown in Figure 6a. It shows that, even as the value of $\beta$ decreases by $10$ times, there is a small change in the number of inhibitors; however, when it increases by $10$ times during a specific period of time, there is a huge difference in the number of inhibitors. Figure 6b demonstrates the impact of the parameter $\epsilon$ on the inhibitors. It shows that, when the value for $\epsilon = 0.43$ decreases by three times, nine times, and so on, the number of inhibitors increases in the different intervals of time and takes more than $120$ days to reduce. And, when the value for the parameter $\epsilon = 0.43$ increases three times, there is a small change in the number of inhibitors; the inhibitors reduced in the same manner as $\epsilon = 0.43$.

Figure 7a illustrates how the nucleators are affected by the parameter $\lambda$. It shows that the number of nucleators increases rapidly and then steadily falls when the value for the parameter $\lambda = 1$ increases by three times. When the value of parameter $\lambda$ decreases by three times, the number of nucleators decreases. In the graph, we can easily see that after a certain interval of time, the number of nucleators decreased in the same way, irrespective of the value of $\lambda$.

The impact of the parameter $\kappa$ on the nucleators is shown in Figure 7b. It demonstrates that when the value of $\kappa = 0.1$ falls by three times, the number of nucleators increases. The graph also demonstrates a significant shift before the peak is reached. Furthermore, while the peak of the nucleators does not significantly change when the value of the parameter is increased by three times, it is slightly displaced.

The effect of the parameter $\upsilon$ on the nucleators can be seen in Figure 8a. It demonstrates that the number of nucleators increases at a specific time interval when the parameter $\upsilon = 17$ falls by three times; it then steadily diminishes in the same way. Additionally, the number of nucleators falls when the value of the parameter $\upsilon$ increases by three and nine times, respectively. Figure 8b demonstrates the impact of the parameter $\eta$ on the nucleators. It shows that as the parameter $\eta = 0.71$ increases by three times, nine times, and more, the number of nucleators decreases but the time interval during which nucleators are present remains unchanged. After a certain time interval, the number of nucleators decreases gradually. Moreover, the number of nucleators increases when the parameter's value is decreased by three times.

Figure 9a shows the change in the parameter $\eta$ affect the rate and degree of mineralization. A three-times decrease in $\eta = 0.71$ resulted in a decrease in the degree of mineralization, whereas a three-fold increase, $9$-fold, and more, resulted in an increase in the degree of mineralization. A change in the value of the parameter $\kappa$ causes a change in the lag time and degree of mineralization, as shown in Figure 9b. A three-fold decrease in $\kappa = 0.10$ resulted in a considerable increment in mineralization lag time; also, a three-fold or more increase merely affected the mineralization lag time. However, a three-fold or more increase in the value of $\kappa$ affects the degree of mineralization. As the value of $\kappa$ increases, the degree of mineralization decreases.

Figure 10a demonstrates that, when $\epsilon = 0.43$ decreases by three-fold, the mineralization lag time decreases and there is a small increase in the degree of mineralization. When the value of $\epsilon$ increases three-fold, the lag time of mineralization increases slowly; however, at $9$-fold, the lag time of mineralization is greater, and there is a decrease in the degree of mineralization. When the value of $\beta = 0.10$ decreases by $10$-fold, as shown in Figure 10b, the lag time decreases and the degree of mineralization increases. A $10$-fold, or more increase results in a greater change in the lag time of mineralization and an approximately 40–50% decrease in the degree of mineralization.

Figure 11a demonstrates that when $\upsilon = 17$ decreases by $3$-fold, there is almost a $2$-fold increase in the degree of mineralization. Furthermore, when $\upsilon$ increases by $3$-fold, the degree of mineralization decreases. There is not much difference in the mineralization lag time due to $\upsilon$. According to Figure 11b, the degree of mineralization decreases when $\lambda = 1$ falls by $3$-fold, leading to greater increases in the degree of mineralization when $\lambda$ increases by 3-fold or more. Due to $\lambda$, there is a minimal difference in the lag time; alternatively, we can say that the difference in mineralization lag time is negligible.

5.   Conclusion

In the current study, a fractional model has been used to investigate the kinetics of bone mineralization, which is a vital bodily function that significantly contributes to the stability of the bones. The local mineral balance and its impact on mineral formation have been predicted dynamically and explained through the defined mathematical model of a nonlinear system of differential equations. Here, a mathematical model was fractionalized by using the C–F fractional derivative with non-singular kernels, where the order of the fractional derivative can be modified to match the experimental data. Numerical simulation was performed using the fractional AB method. It has been ascertained that the presence of minerals can be depicted at an early stage with the help of a fractional model which is not true for an integer-order model. Sensitivity analysis has been performed to study the effects of various parameters, i.e., the collagen cross-linking rate $\kappa$, rate of release of inhibitors $\beta$, rate of removal of inhibitors $\epsilon$, number of nucleators per mature collagen molecule $\lambda$, rate of removal of nucleators $\upsilon$, and rate of mineralization $\eta$. In particular, a significant impact has been observed on the mineral concentration level in relation to a variety of parameters. The capacity of the suggested fractional model to forecast the mineralization at an early stage can help researchers and decision-makers to take preventive actions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Researchers would like to thank the Deanship of Scientific Research, Qassim University for funding publication of this project.

Conflict of interest

The authors declare there is no conflict of interest.