*Article* **Artificial Neural Network Solution for a Fractional-Order Human Skull Model Using a Hybrid Cuckoo Search Algorithm**

**Waseem1,\*, Sabir Ali 2, Shahzad Khattak 2,3, Asad Ullah 4,5, Muhammad Ayaz 6, Fuad A. Awwad <sup>7</sup> and Emad A. A. Ismail <sup>7</sup>**


**Abstract:** In this study, a new fractional-order model for human skull heat conduction is tackled by using a neural network, and the results were further modified by using the hybrid cuckoo search algorithm. In order to understand the temperature distribution, we introduced memory effects into our model by using fractional time derivatives. The objective function was constructed in such a way that the *L*2−error remained at a minimum. The fractional order equation was then calculated by using the proposed biogeography-based hybrid cuckoo search (BHCS) algorithm to approximate the solution. When compared to earlier simulations based on integer-order models, this method enabled us to examine the fractional-order (FO) cases, as well as the integer order. The results are presented in the form of figures and tables for the different case studies. The results obtained for the various parameters were validated numerically against the available literature, where our proposed methodology showed better performance when compared to the least squares method (LSM).

**Keywords:** boundary value problems; fractional derivatives; heat conduction; BHCS algorithm; Cuckoo search; numerical method; human head

#### **1. Introduction**

The use of electronic devices is increasing day by day. One reason behind this is the advancement of technology and its applications in various sectors. This advancement in technological equipment has some side effects, especially when it crosses some limit in its use. Some of the devices and systems are Bluetooth, mobile phones, and other headphone-like devices. These devices produce thermal waves that pass through the skin and enter the human head, damaging various tissues, including the brain. The use of these electronic devices produces brain damage and other neural disorders, as explained in the references [1–3]. Furthermore, the analysis of the effects of the thermal and nonthermal waves are numerically and experimentally analyzed by Bernardi et al. [4]. The flow symmetry of heat is widely analyzed by many researchers due to its experimental and theoretical applications [5]. The energy transfer from the electronic object to the human head follows the one-sided symmetry in the skull. The facial, brain, and skull symmetries are well explained by Ratajczak et al. [6].

**Citation:** Waseem; Ali, S.; Khattak, S.; Ullah, A.; Ayaz, M.; Awwad, F.A.; Ismail, E.A.A. Artificial Neural Network Solution for a Fractional-Order Human Skull Model Using a Hybrid Cuckoo Search Algorithm. *Symmetry* **2023**, *15*, 1722. https://doi.org/10.3390/ sym15091722

Academic Editors: Francisco Martínez González, Mohammed K. A. Kaabar and Junesang Choi

Received: 12 July 2023 Revised: 26 August 2023 Accepted: 5 September 2023 Published: 8 September 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The analysis of heat in a human head became famous after the work of Flesch [7]. In this work, Flesch used a differential equation for the analysis of heat. This work was further examined by Gray [8] in 1980, which provides a theoretical approach to heat transfer analysis in terms of the human head. The human skull produces more heat on the outer layer than the center. On the other hand, when the surrounding temperature is reduced, the heat production is more peripheral to the skull. Simply put, the temperature outside and the radial distance from the skull's center affect how much heat is produced inside the human skull. Anderson and Arthurs [9] used the complementary extremum approach to analyze this famous problem. Makinde [10] analyzed the human skull problem by using the non-perturbation approach. Raja et al. [11] implemented the stochastic approach to study the human skull problem. Abdelhakem and Youssri [12] analyzed the Lane-Emden and Bratu equations by using the spectral Legendre's algorithm. Youssri et al. [13] used the wavelets approach for a solution of the Lane–Emden equations. A more brief analysis using the numerical approaches to the solution of D.Eqs. can be found in the literature [14,15].

The methodology and the model modifications are both points of interest to researchers. In recent years, the use of fractional derivatives in differential equations has been widely implemented [16,17]. The applications of fractional derivatives are well explained by Podlubny [18]. The use of fractional derivatives in the field of differential equations is explained by Aleksandrovich et al. [19]. Wang [20] studied the febrifuge effect for analyzing the fractional-order human skull problem. The concept of Caputo-type derivatives was introduced by Kumar et al. Kumar et al. [21] for use in differential equations. At the same time, the concept of Caputo-type derivatives for delay-type differential equations was introduced by Odibat et al. [22]. The concept of fractional derivatives in applications can be found in ecology [23,24], psychology [25], chemistry [26], epidemiology [27–30], and physics [31]. Yavuz and Sene [32] examined how various parameters affect fractional-order second-grade fluid flow. Hammouch et al. [33] numerically simulated a fractional chaotic system with changing order. Yavuz [34] studied the classic and generalized Mittag-Leffler kernels used in the fractional derivative definition in the European option pricing model.

The applications of fractional derivatives are not limited to a single definition. The solution to Cauchy and Dirichlet problems are studied by Avci et al. [35] by using the Caputo-Fabrizio derivative definition. Erturk et al. [36] developed a unique Caputo fractional derivative for the corneal shape model of the human eye. The recent literature shows the application of fractional calculus, which can be found in the following references [37–39].

The following is the integer-order model for temperature distribution in the human skull:

$$\begin{array}{l} T''(r) + \frac{2T'(r)}{r} + \lambda.exp(-mT) = 0, \\ T'(0) = 0, \ T'(1) = \mathcal{N}\_B(1 - T). \end{array} \tag{1}$$

Here, *T*, *r*, *NB*, *m*, and *λ* denote the temperature, radial distance, Biot number, metabolic thermogenesis slope, and thermogenesis heat production, respectively.

We introduce the Riemann–Liouville definition of the fractional-order derivative of the function ∈ *<sup>C</sup><sup>d</sup>* <sup>−</sup><sup>1</sup> below [18]:

$$D\_t^{\mu} \mathbb{S}(t) = \begin{cases} \frac{d^{\gamma} \mathbb{S}(t)}{d\mathbb{S}^{\gamma}}, & \mu = \gamma \in \mathbb{Z},\\ \frac{1}{\gamma(\gamma - \mu)} \int\_0^t (t - \sigma)^{\gamma - \mu - 1} \mathbb{S}^{\gamma}(\sigma) d\sigma, & \gamma - 1 < \mu < \gamma, \gamma \in \mathbb{Z}. \end{cases} \tag{2}$$

In light of Equation (2), the suggested classical BVP (1) is transformed into a fractionalorder generalized form:

$$\begin{array}{ll} \frac{2}{r} \, ^c D\_0^\nu T(r) + ^c D\_0^\mu T(r) + \lambda .exp(-mT) = 0, \\ \qquad T'(0) = 0, \; \, T'(1) = \mathcal{N}\_B(1 - T). \end{array} \tag{3}$$

Here, *r* ∈ [0.1, 1] and the derivatives of the function *T*(*r*) with fractional orders <sup>0</sup> <sup>&</sup>lt; *<sup>ν</sup>* <sup>≤</sup> 1 and 1 <sup>&</sup>lt; *<sup>μ</sup>* <sup>≤</sup> 2, are represented by the symbols *cD<sup>ν</sup>* <sup>0</sup> and *cD<sup>μ</sup>* <sup>0</sup> respectively.

The graphical abstract of the proposed methodology is given in Figure 1, whereas the structure of the neural network is presented in Figure 2.

**Figure 1.** Graphical abstract of the given model.

**Figure 2.** Graphical presentation of the ANN structure for the given model.

#### *Our Contribution*

The primary objective of this research study is to develop an approximate solution for a fractional-order human head heat conduction model by using the BHCS algorithm. More specifically, we summarize our contributions as follows:


In this article, the proposed methodology of BHCS is discussed in Section 2, and the numerical step-up and various proposed cases (with graphs and tables) are presented in Section 3 and are discussed in detail. At the end, a conclusion is provided in Section 4 of the article.

#### **2. The Proposed Methodology**

The approximate solution (*T*ˆ(*r*)) of the fractional-order model for temperature in the human head by using feed-forward neural networks with the help of an exponential function is given as

$$\hat{T}(r) = \sum\_{i=1}^{m} \alpha\_i e^{\omega\_i r + \beta\_i} \,\,\,\,\,\,\tag{4}$$

$$\frac{d^\mu \hat{T}(r)}{dr^\mu} = \sum\_{i=1}^m \alpha\_i r^{-\mu} e^{\beta\_i} E\_{1,1-\mu}(\omega\_i r),\tag{5}$$

$$\frac{d^\nu \hat{T}(r)}{dr^\nu} = \sum\_{i=1}^m \alpha\_i r^{-\nu} e^{\theta\_i} E\_{1,1-\nu}(\omega\_i r) \tag{6}$$

where *αi*, *ω<sup>i</sup>* and *β<sup>i</sup>* are the weights given as *α* = [*α*1, *α*2, ... , *αm*], *ω* = [*ω*1, *ω*2, ... , *ωm*], *<sup>β</sup>* = [*β*1, *<sup>β</sup>*2, ... , *<sup>β</sup>m*], and *<sup>m</sup>* represents the number of neurons. Moreover, *<sup>d</sup>νT*ˆ(*r*) *dr<sup>ν</sup>* and *<sup>d</sup>μT*ˆ(*r*) *dr<sup>μ</sup>* are the fractional derivatives of the series solution.

#### *2.1. Fitness Function*

In the fitness function, we compute the absolute error and make an optimization process to minimize the error , *i.e.,* when <sup>→</sup> 0, then *<sup>T</sup>*ˆ(*r*) <sup>→</sup> *<sup>T</sup>*(*r*).

The fitness function for the transformed equations is given as

$$
\mathfrak{e} = \mathfrak{e}\_1 + \mathfrak{e}\_2.\tag{7}
$$

where <sup>1</sup> represents the mean squared error for a given differential equation (DE), and <sup>2</sup> represents the conditions on it. Therefore, we have

$$\varepsilon\_{1} = \frac{1}{k} \sum\_{k=1}^{K} \left( \frac{2}{r\_{k}} \, ^{c}D\_{0}^{\nu} \hat{T}\_{k} + ^{c}D\_{0}^{\mu} \hat{T}\_{k} + \lambda.exp(-m\hat{T}\_{k}) \right)^{2},\tag{8}$$

and

$$\varepsilon\_{2} = \frac{1}{2} \left( (\hat{T}\_{0}')^{2} + (\hat{T}\_{1}' - \text{N}\_{B}(1 - \hat{T}))^{2} \right), \tag{9}$$

where *h* = *<sup>u</sup> <sup>k</sup>* , *<sup>T</sup>*<sup>ˆ</sup> *<sup>k</sup>* = *T*ˆ(*rk*), and *rk* = *kh*.

#### *2.2. Cuckoo Search (CS) Technique*

The cuckoo search (CS) algorithm follows the breeding behavior of the cuckoo bird [41]. In this algorithm, other birds give their eggs to others' nearest nests. When the host bird finds it, she adopts two methods: either to remove the eggs or to find a new nearest nest to lay their own eggs. In this process, the host bird's eggs indicate a solution, and the cuckoo bird's eggs display a fresh potential resolution [42].

This procedure is explained in [43]:


We assume, *yi* = *yi*1, *yi*2, *yi*3, ... , *yiD* as *i th* egg positions. The egg is defined as a solution, and Lévy flights update the new solution *ynew <sup>i</sup>* as follows:

$$y\_i^{new} = y\_i^{old} + a(y\_l - x\_\mathcal{g}) \oplus \text{Levy}(\mathcal{g}),\tag{10}$$

$$y\_i^{new} = y\_i^{old} + \frac{0.01u}{|v|^{\frac{1}{\beta}}}(y\_i - y\_{\mathcal{S}})\_\prime \tag{11}$$

where ⊕ is the entry-wise product, *β* indicates the Lévy flight exponent, *α* > 0 is the cuckoo's step size, *xg* is the optimal sample, and *u* and *v* are random numbers. Furthermore,

$$v \sim N(0, \sigma\_v^2), \qquad u \sim N(0, \sigma\_u^2), \tag{12}$$

$$\sigma\_{\rm u} = \left[ \frac{\sin \frac{\pi \beta}{2} \cdot \Gamma(1+\beta)}{2^{\frac{\beta-1}{2}} \beta \cdot \Gamma(\frac{1+\beta}{2})} \right]^{\frac{1}{\beta}}, \qquad \sigma\_{\rm v} = 1. \tag{13}$$

Here, the function *σ<sup>u</sup>* is controlled by the parameter *β* and the Γ function. In CS, the discovered nests are replaced using a discovery operator that takes the probability *pa* into account. Thus, we have the updated solution given as follows:

$$y\_{ij}^{\text{new}} = \begin{cases} y\_{ij}^{\text{old}} + \text{rand} \cdot (y\_{r1}j(k) - y\_{r2}j(k)) & \text{if } P > p\_{a\nu} \\ y\_{ij}^{\text{old}}(k) & \text{else}. \end{cases} \tag{14}$$

Here, *y*new *ij* represents the *<sup>j</sup>*th component of the *<sup>i</sup>*th solution *<sup>y</sup>*new *<sup>i</sup>* , *yr*1,*<sup>j</sup>* is the *j*th element of the solution *yr*1, and *yr*2,*<sup>j</sup>* is *j*th element of the solution *yr*2. Moreover, *r*1 and *r*2 are two distinct integers within in [1, *NP*], where *NP* denotes the size of the population, and *pa* is the discovery denoting probability.

#### *2.3. Biogeography-Based Optimization*

An evolutionary biogeography-based optimization method (BBO) was motivated by several traits of animals found on islands. In BBO, NP habitats (solutions) are used to randomly initialize the population. Each generation ranks the population from best to worst.

Here, we define *λ*ˆ and *μ*ˆ as the particular habitat's immigration and emigration rates, as given in [44]:

$$\begin{cases} \hat{\lambda}\_i = I \left( 1 - \frac{\mathcal{S}\_i}{N\mathcal{P}} \right), \\ \hat{\mu}\_i = E \frac{\mathcal{S}\_i}{N\mathcal{P}}. \end{cases} \tag{15}$$

where *E* = *I* = 1 are the immigration rates, *s*ˆ*<sup>i</sup>* is a species of a certain population, which is defined as {*s*ˆ*<sup>i</sup>* = *NP* − *i*}, *i* ∈ *N*. The changing parameter updates the corresponding solution, and BBO also utilizes the mutation operator to update the solution accordingly.

#### *2.4. Hybrid Cuckoo Search*

In order to further improve the best nests obtained from the CS, we applied BBO. Both exploration and exploitation were employed alternatively. By combining exploration and exploitation, the BBO-based heterogeneous cuckoo search (BHCS) method was designed as a hybrid meta-heuristic. The proposed BHCS algorithm comprises two primary steps: heterogeneous CS and biogeography-based discovery.

#### The Methodology of Heterogeneous CS

The first component of BHCS employs the Lévy flights and a quantum mechanismbased heterogeneous CS. Heterogeneous CS is based on quantum mechanics [42,45].

$$y\_i^{\text{new}} = \begin{cases} y\_i^{\text{old}} + \mathfrak{a} \cdot \left( y\_i - y\_\mathcal{g} \right) \oplus \text{Lévy}(\beta) & \frac{2}{3} < s\_r < 1, \\ \bar{y} + L \cdot \left( \bar{y} - y\_i^{\text{old}} \right) & \frac{1}{3} < s\_r \le \frac{2}{3}, \\ y\_i^{\text{old}} + \varepsilon \cdot \left( y\_\mathcal{g} - y\_i^{\text{old}} \right) & \text{else}. \end{cases} \tag{16}$$

Here, the terms *L* = ln ( <sup>1</sup> *<sup>η</sup>* ), = *<sup>δ</sup>eη*, and *yg* refer to the iteration's best solution, *sr* is the number *<sup>η</sup>* <sup>∈</sup> [0, 1], and *<sup>y</sup>*¯ <sup>=</sup> <sup>1</sup> *NP* <sup>∑</sup>*NP <sup>i</sup>*=<sup>1</sup> *yi* is the average of the solutions. Equation (16) demonstrates that three equations are used in a heterogeneous cuckoo search to update the answers with the exact probabilities. The 1st equation is derived from the Lévy flights in the original CS, and the 2nd and 3rd equations are used to update the results by using the quantum-based algorithm. The search space is diversified by updating the solutions with heterogeneous rules, which move toward the actual global region.

#### **3. Results and Discussion**

In this section, we calculate four individual case studies and compute the results using our proposed BHCS as a global search technique. The case studies are shown in Table 1. A total of 100 independent runs were performed for each case study by taking the domain *r* ∈ [0.1, 1] with 0.1 step size.

**Table 1.** Different cases with the variation of fractional derivative parameters.


The formulation for these case studies is given below:

#### *3.1. Case Study 1*

For this case, the fitness function with the boundary conditions is given as

$$\varepsilon\_{1} = \frac{1}{11} \sum\_{k=1}^{11} \left( \frac{2}{r\_{k}} .^{c}D\_{0}^{1}\hat{T}\_{k} + ^{c}D\_{0}^{2}\hat{T}\_{k} + \lambda.exp(-m\hat{T}\_{k}) \right)^{2} . \tag{17}$$

$$\varepsilon\_2 = \frac{1}{2} \left( (\hat{T}\_0')^2 + (\hat{T}\_1' - \mathcal{N}\_B(1 - \hat{T}))^2 \right). \tag{18}$$

We considered the parameters *λ* = 1, *m* = 1, *Nb* = 1. By using BHCS, the best weights for this case are given in the following equation.

$$\mathcal{T}\_{\mathbf{c}\_{1}} = \begin{cases} 1.4025e^{( -0.2030r - 0.7130)} - 0.0856e^{( -1.4281r + 0.5922)} \\ \quad + 0.2617e^{( -0.0189r - 0.1831)} - 0.7062e^{( 0.5749r - 1.612)} \\ \quad - 2.4658e^{( -0.5143r - 1.7698)} + 1.7948e^{( -0.0388r - 0.7866)} \\ \quad - 0.1858e^{( -0.5612r - 0.5207)} + 0.7182e^{( -1.3911r - 1.2173)} \\ \quad + 0.9790e^{( 0.3547r - 2.6460)} - 0.0509e^{( -1.6107r - 1.0547)} \end{cases} \tag{19}$$

Here, Equation (19) is a series solution for case study 1.

#### *3.2. Case Study 2*

The fitness function for the current case study is formulated below:

$$\varepsilon\_{1} = \frac{1}{11} \sum\_{k=1}^{11} \left( \frac{2}{r\_{k}} .^{c}D\_{0}^{0.70} \hat{T}\_{k} + ^{c}D\_{0}^{1.70} \hat{T}\_{k} + \lambda.exp(-m\hat{T}\_{k}) \right)^{2} . \tag{20}$$

$$\varepsilon\_2 = \frac{1}{2} \left( (\hat{T}\_0')^2 + (\hat{T}\_1' - N\_B(1 - \hat{T}))^2 \right). \tag{21}$$

Here, we take the parameters *λ* = 1, *m* = 1, *Nb* = 1. The best weights for this case are given as

$$\hat{T}\_{c\_2} = \begin{cases} -0.7538e^{( -0.8125r + 0.7228)} + 1.3636e^{( -0.3377r - 0.4151)} \\ + 0.8281e^{( -0.2094r - 1.3968)} + 0.4684e^{( -0.4236r - 0.2732)} \\ + 0.9507e^{( -1.3841r - 1.4493)} + 1.1338e^{( -1.1852r - 0.9906)} \\ -0.6161e^{( -5.0000r - 4.1270)} - 0.1282e^{( -0.7593r - 1.3885)} \\ + 1.9087e^{( -0.3513r - 0.8274)} + 0.0553e^{( 0.1238r - 0.1236)} \end{cases} \tag{22}$$

Here, the above equation is the corresponding series solution for this special case 2. The absolute errors (AE) are presented in Table 2. The approximate solutions of the given fractional model, which were obtained by using the series solution of Equation (19), are illustrated in Figure 3. The best fitness values are evaluated by considering the above conditions, and the results are displayed in Figure 4.


**Table 2.** Minimum Absolute Errors ().

We analyzed the FO heat conduction model for the human head given in Equation (3) by using an interval of [0.1, 1]. Four different cases were considered by choosing varying values of *μ*, *ν* for the fixed parameters *λ* = 1, *Nb* = 1. The approximate solutions for both cases, 1 and 2, are presented in Figure 3. These graphs show that when we decrease the order of the fractional parameters from an integer order to a non-integer, the temperature profile jumps to 1.17 from 1.16. A decreasing trend is observed when *r* −→ 1. This trend is more sharp in the integer order when compared to the fractional order. The fitness of the functions is shown in Figure 4, where the red lines show the mean. Almost all the values are bounded by the box, and the distance from the mean positions is displayed on the vertical line.

**Figure 3.** Graphical representation of the solutions for (**a**) Case study 1 and (**b**) Case study 2.

**Figure 4.** Graphical representation of fitness functions for (**a**) Case study 1 and (**b**) Case study 2.

#### *3.3. Case Study 3*

In case 3, we considered *μ* = 1.80 and *ν* = 0.80 by choosing *λ* = 1, *m* = 1, *Nb* = 1. The fitness functions for this case are given by

$$\epsilon\_1 = \frac{1}{11} \sum\_{k=1}^{11} \left( \frac{2}{r\_k} .^{\circ}D\_0^{0.80} \hat{T}\_k + ^{\circ}D\_0^{1.80} \hat{T}\_k + \lambda.exp(-m\hat{T}\_k) \right)^2,\tag{23}$$

$$\epsilon\_2 = \frac{1}{2} \left( (\hat{T}\_0')^2 + (\hat{T}\_1' - N\_B(1 - \hat{T}))^2 \right). \tag{24}$$

The corresponding best weights are given below:

$$\hat{T}\_{\varepsilon\_3} = \begin{cases} 0.3527e^{(0.7055r - 2.1888)} - 1.1520e^{(-0.7899r - 1.8388)} \\ \quad + 1.4121e^{(-3.4830r - 1.8823)} + 1.3158e^{(-0.4966r - 0.9646)} \\ \quad + 0.1036e^{(0.3872r - 0.0793)} - 1.0951e^{(-0.9004r - 0.8300)} \\ \quad - 0.2317e^{(-3.6303r - 0.1507)} + 0.1977e^{(-1.2789r - 0.2822)} \\ \quad + 0.8166e^{(-0.0468r + 0.4276)} - 0.5767e^{(0.5119r - 0.9304)} . \end{cases} \tag{25}$$

#### *3.4. Case Study 4*

For this case study, we took *μ* = 1.90 and *ν* = 0.90. So, the corresponding fitness functions take the following form:

$$\varepsilon\_{1} = \frac{1}{11} \sum\_{k=1}^{11} \left( \frac{2}{r\_{k}} \prescript{c}{}{D\_{0}^{0.90}} \hat{\mathbf{T}}\_{k} + \prescript{c}{}{D\_{0}^{1.90}} \hat{\mathbf{T}}\_{k} + \lambda . exp(-m\hat{\mathbf{T}}\_{k}) \right)^{2} \tag{26}$$

$$\epsilon\_2 = \frac{1}{2}\left( (\hat{T}\_0')^2 + (\hat{T}\_1' - N\_B(1 - \hat{T}))^2 \right). \tag{27}$$

By choosing *λ* = 1, *m* = 1, *Nb* = 1, we have

$$
\hat{T}\_{c4} = \begin{cases}
+ 0.1558e^{(-1.3456r - 1.2083)} - 0.9332e^{(0.3765r - 0.7752)} \\
+ 1.3882e^{(-0.1069r + 0.1492)} - 0.3297e^{(-0.4075r - 1.2062)} \\
\end{cases} \tag{28}
$$

Here, *T*ˆ *<sup>c</sup>*<sup>3</sup> and *T*ˆ *<sup>c</sup>*<sup>4</sup> are the series solutions for cases 3 and 4. The AE are plotted in Table 2, whereas the approximate solutions are presented in Figure 5. The fitness values are given in Figure 6. The values of AE for all the case studies are presented in Table 2.

Similarly, in cases 3 and 4, when we increase the fractional parameters that nearly approach the integer, the solution profiles fall from 1.65 to 1.61. As a result, the suggested fractional-order graph, which takes into account the radial distance (*r*), Biot number (*NB*), metabolic thermogenesis slope parameter (*m*), and thermogenesis heat production parameter (*λ*), provides a more accurate representation of the distribution of temperature within the human skull.

The fitness functions for cases 3 and 4 are displayed in Figure 6. The horizontal red line shows the mean, and the red addition symbols show the positions of the results from this point. In both cases, the results are in the range of 10−<sup>5</sup> and 10−4, respectively. This further recommends that the fitness functions remain as minimal as possible.

In Table 2, the results for the minimum values of the absolute error are displayed numerically. These cases are chosen in such a way that the deviations from the integer order to a fractional order are clearly visible as time varies. First, the decreasing trend from the integer order is observed for various fractional parameters. In the last two columns, the increasing trend towards the integer order is displayed. If we compare the results of the second case and the fourth, we see that the results initially go toward the worst and

then beat the integer order at *r* = 0.7. Similarly, the BHCS results are compared against the available literature in Table 3. In case 1, the results for the integer order are almost the same as the LSM. In cases 2, 3, and 4, the results of BHCS are nearer the integer-order solution when compared to LSM. This proves the validity and efficiency of our proposed method, BHCS.

**Table 3.** Comparison of the approximate solutions of the BHCS neural network using the least squares method (LSM) [40].


**Figure 5.** Graphical representation of the solutions for (**a**) Case study 3 and (**b**) Case study 4.

**Figure 6.** Graphical representation of fitness function for (**a**) Case study 3 and (**b**) Case study 4.

#### **4. Conclusions**

In this article, we discussed the distribution of temperature within the human skull by considering the fractional derivative. In order to solve the proposed model, we utilized the biogeography-based hybrid cuckoo search (BHCS) algorithm to then be used on the transformed fractional-order equation. The following are the main features obtained based on our analysis.


**Author Contributions:** Conceptualization, W. and S.A.; methodology, S.K.; software, W. and A.U.; validation, M.A., E.A.A.I., S.K., F.A.A. and A.U.; formal analysis, W. and E.A.A.I.; investigation, S.A.; resources, A.U.; data curation, W.; writing—original draft preparation, M.A.; visualization, A.U.; supervision, W.; review writing, E.A.A.I., F.A.A., W., S.A. and A.U.; funding, E.A.A.I. and F.A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was funded by King Saud University, Riyadh, Saudi Arabia.

**Data Availability Statement:** No data are used in this study.

**Acknowledgments:** Researchers' supporting project number: (RSPD2023R1060), King Saud University, Riyadh, Saudi Arabia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
