Construction a distributed order smoking model and its nonstandard finite difference discretization

1.   Introduction

In recent years, fractional calculus has become a popular field of study because of its effective application in different scientific fields, such as statistics, applied mathematics, dynamics, mathematical biology, control theory, optimization, and chaos theory. Using distributed order differential equations, which is a general form of both ordinary and fractional-order differential equations, is very important for mathematical modeling. Also, distributed order differential equations are important in biology and engineering models, as they include the density function. The selection of the density function allows determining the dynamics of the models used in different situations.

Caputo introduced the idea of distributed order fractional derivative and informed about the usage areas of differential equations with this type of derivative ^[8]. Caputo used the definition of distributed order derivative in filter transfer functions and physical applications ^[9]. Then Hartley and Lorenzo used continuously distributed order equations for fractional systems ^[17]. Bagley and Torvik discussed the existence and solution of distributed order differential equations ^[5,6]. Caputo introduced distributed order differential equations in modeling dielectric induction and diffusion ^[10].

Recently, the studies of distributed order differential equations have attracted great attention from many researchers. Ford and Morgado studied the existence and uniqueness of solutions of distributed order differential equations ^[16]. Luchko proved the uniqueness and continuous dependence of the distributed sequential time-fractional diffusion equation based on the initial conditions ^[25].

Refahi et al. focused on analytical solutions of distributed order fractional linear equation systems in ^[35]. Katsikadelis introduced efficient numerical methods for distributed order fractional differential equations in ^[20]. Diethelm and Ford presented the solutions to distributed order differential equations numerically, including error analysis ^[12]. Li and Wu presented a numerical method for the distributed order diffusion equation in ^[22]. Li and Wu used the classical quadrature formula, multi-term fractional diffusion equations, and the kernel method. Morgado and Rebelo developed an implicit scheme for a numerical approximation of the distributed order reaction-diffusion equation ^[31].

Studies on stability analysis of distributed order differential equations were presented by Aminikhah et al. and Najafi et al. ^[3,4,37]. Amin et al. studied the distributed order time-fractional equations solutions using the Haar wavelet method ^[2].

In this study, we aimed to obtain a numerical solution of a new distributed order fractional smoking model. Because, by choosing the density function, we can find interpretations for different situations in the smoking model. Hence, the effect of the smoking model on different parameters can be examined with these comments.

Smoking affects almost all organs in our body and causes diseases. Smoking is harmful to the lungs, and it destroys the vesicles (alveoli), the smallest unit where oxygen enters the lungs. Therefore, it causes diseases such as pneumonia, asthma, and tuberculosis. Smoking also causes significant damage to the immune system. So, it is very difficult for a smoker to overcome diseases such as Coronavirus Disease (COVID), Severe Acute Respiratory Syndrome (SARS), and Middle East Respiratory Syndrome (MERS). That's why it is important to obtain data for smoking dynamics in different situations. For these reasons, mathematical models related to the effects of smoking have gained importance, especially after 2000.

First, Castillo-Garsow et al. proposed a simple mathematical model for smoking cessation ^[11]. They envisioned a system with a total fixed population divided into three classes: potential smokers, former smokers, and smokers. Following this study, Sharomi and Gumel developed a smoking model by introducing the light and chain classes ^[38]. Zaman improved the work of Castillo-Garsow and developed a model that took into account the occasional smokers quotient in the smoking cessation model ^[44]. After these studies, many articles have been written on stability analysis, mathematical analysis, and approximate solutions on the smoking model. For detailed information, you can review ^{[15,18,19,24,40,41,42]} resources.

The rest of this paper is organized as follows. In Section 2, we present some definitions of distributed order differential equations. In Section 3, we define and discretize a new distributed order smoking model. Also here, some important information gives about the NSFD method and stability analysis of the system. In Section 4, we show numerical simulations graphically. Finally, in Section 5, we discuss the results and conclude the paper.

2.   Basic definitions about distributed order differential equations

In this section, we give basic information and definitions about distributed order and fractional-order differential equations. Then, information is given about stability analysis for distributed order fractional equation systems. These are the important books and papers related to this topic ^{[7,8,21,23,34,37]}.

Definition 2.1. ^[34] Riemann Liouville fractional derivatives of order $\alpha$ is defined by

where $f(t)$ can be integrated in the $[a, b]$, $n-1 < \alpha \leq n \ \ (n \in \mathbb{N^+})$ and $\Gamma()$ is a Gamma function.

Definition 2.2. ^[34] As in Definition (2.1), let $f(t)$ can be integrated in the $[a, b]$ and $n-1 < \alpha \leq n \ \ (n \in \mathbb{N^+})$, Caputo fractional derivative is defined by

Definition 2.3. ^[34] Suppose that the function $f ^ {(z)} (t)$ is continuous in the closed range $[a, t]$ and has derivatives for $z = 1, 2, 3, ..., n + 1$. Let $n$ be an integer and $p$ satisfy $n < p < n + 1$. In this case, Grünwald-Letnikov fractional derivative definition is given as:

or

Definition 2.4. ^[8] Caputo first defined an integral operator for distributed order differential equations. This operator is defined by

where $\alpha \in (0, 1]$ and $\int _{0} ^{1} w(\alpha) = k > 0$. Additionally, $D^{\alpha}_{t} f(t)$ is a fractional derivative function and it can be chosen as Riemann-Liouville, Caputo or Grünwald-Letnikov ^[8,9,10]. Another important function $w (\alpha)$ is the density function of the distributed order differential equation. The reason why these distributed equations are considered the general form of ordinary and fractional differential equations is the importance of density function selection.

Distributed order differential equation is defined by

Definition 2.5. ^[26] The approximate Grünwald-Letnikov formula for distributed order differential equations is defined by

By rearranging in (2.7), this formula can be expressed as:

where $q_i ^{\alpha} = (1-\frac{1+\alpha}{i}) q_{i-1} ^{\alpha}$, $q_0 ^{\alpha} = h^{-\alpha}, i = 0, 1, 2, 3, ..., n$ and $h$ is assumed to be very small ^[14,43].

2.1.   Some definitions about stability analysis

Let

where $u(t) = (x_1(t), x_2(t), ..., x_n(t))^T$ and the function $G$ is differentiable.

A general numerical scheme with a step size $h$, that approximates the solution $u(t_k)$ of the system (2.9) can be defined in the form:

where $D_h (u_k) \approx (\frac{dx_1}{dt}, \frac{dx_2}{dt}, ..., \frac{dx_n}{dt})^T$, $u_k \approx u(t_k)$, $t_k = t_0+kh$ and $F_h(G; u_k)$ approximates the right-hand side of system (2.9).

Definition 2.6. ^[13] Let $E^*$ an equilibrium point and $J(E^*)$ the Jacobian of system (2.9). An equilibrium point $E^*$ is called linear stable if $Re(\lambda) < 0$ for all $\lambda \in \sigma(J(E^*))$ and linear unstable if $Re(\lambda) > 0$ for at least one $\lambda \in \sigma(J(E^*))$.

Lemma 2.7. ^[13]Assume that (2.10) has the following form:

where the function $H$ is differentiable. A fixed point $E^*$ is stable if and only if all eigenvalues of $J(E^*)$ are less than one in absolute values.

While analyzing the eigenvalues of $J(E^*)$, we will use Schur-Cohn Test. The generalized Schur-Cohn test can be expressed in the following form. This test is an algebraic criterion to decide whether a discrete system is stable or not ^[36]. Characteristic equation can be expressed as:

where $a_n > 0$. Some of the determinants to be used are expressed as follows:

The Schur-Cohn test states that the equilibrium point is stable only if all the following conditions are satisfied:

$1) \ P(1) > 0,$

$2) \ (-1)^n P(-1) > 0,$

$3) \ |a_0| < a_n,$

$4) \ |c_{i, 0}| > |c_{i, n+1-i}|, \ \ \ i = 2, 3, ..., n-1.$ ^[36]

The conditions obtained by the Schur-Cohn test are given below.

Lemma 2.8. ^[36]For the quadratic equation $\lambda^2+a_1 \lambda +a_0 = 0, \ \ (a_2 = 1)$, all roots satisfy $|\lambda_i| < 1, \ \ i = 1, 2$ if and only if the following three conditions are satisfied :

(i) $1+a_1+a_0 > 0$,

(ii) $1-a_1+a_0 > 0$,

(iii) $a_0 < 1$.

Since the value of n for the smoking model is 5, the lemma to be used in this case is as follows:

Lemma 2.9. ^[36] If $P(\lambda) = \lambda^5+a_{4} \lambda^{4}+a_{3} \lambda^{3}+a_{2} \lambda^{2}+a_1 \lambda +a_0 = 0\ \ (a_5 = 1)$, all roots satisfy $|\lambda_i| < 1, \ \ i = 1, 2, 3, 4, 5$ if and only if the following six conditions are satisfied:

(i) $1+a_4+a_3+a_2+a_1+a_0 > 0$,

(ii) $1-a_4+a_3-a_2+a_1-a_0 > 0$,

(iii) $a_0 < 1$,

(iv) $|c_{2, 0}| > |c_{2, 4}|$,

(v) $|c_{3, 0}| > |c_{3, 3}|$,

(vi) $|c_{4, 0}| > |c_{4, 2}|$.

3.   A novel distributed order smoking model and its discretization

The relationship between simple mathematical modeling and biological or physical system, integer-order differential equations show the dynamics of systems. Integer-order differential equations conjugate the relationship between complex system parameters in mathematical modeling and describe the variation of structure within them, nonlinearity, and multi-scale behavior.

Fractional calculus has attracted a significant amount of attention by researchers in recent years and different aspects of the subject have been investigated. This is because the fractional derivative is an important tool for explaining the dynamic behavior of various physical systems. The strength of these differential operators is their nonlocal property not found in integer- order differential operators. The distinguishing features of fractional differential equations are that they summarize the memory and transmitted properties of many mathematical models. Fractional-order models are more realistic and practical than conventional integer-order models. Arbitrary order derivatives are powerful tools for the appreciation of the dynamic behavior of various biomaterials and systems. The most recurring feature of these models is their global feature, which is not present in classical layout models ^[18].

Since the smoking model is also a model of the above-mentioned styles, it is important to better determine the dynamics and reach more realistic values. For this reason, we have used distributed order differential equations in this article.

In this section, the ordinary and fractional smoking model defined by Singh et al. in ^[40] is generalized. This generalization is carried out by using the distributed order differential equation. The new system for the smoking model is defined as follows:

In this system, the whole population is divided into five different subgroups. These subgroups mean P(t): Potential smokers, L(t): Occasional smokers, S(t): Heavy smokers, Q(t): Temporary quitters and R(t): Smokers who quit permanently. The expressions represented by the meanings of the constants are given in Table 1.

$\bf{Reproductive \ \ number:}$ Finding boundary conditions to control the status of the population is very important in such endemic systems. $R_o$ value, also called basic reproduction number, is needed in these cases. Two different situations can be expressed according to the $R_o$ value.

i) In smoking-free equilibrium point case, if $R_o < 1$, which indicates that the disease will die.

ii) In smoking present equilibrium point case, if $R_o > 1$, it means the disease will spread.

To find the value of $R_o$, the Jacobian $(J)$ matrix must be expressed in the type $J = U-V$. Then matrix $Z = UV^{-1}$ must be obtained. Using this matrix $Z$, the determinant $|Z-\lambda I|$ is found and the $\lambda$ value gives us the $R_o$ value. In this smoking model, $R_o = \frac{\tau_1 (1-\tau_2) \gamma}{(\rho+\tau_1)(\rho+\gamma)}$ ^[1]. For detailed information on the importance and analysis of the $R_o$ value in endemic models, see Shi et al. ^[39].

Numerical methods such as Runge-Kutta, Adams methods, and Theta methods based on finite difference approaches are often used to study the dynamics of interacting populations. However, the disadvantages of these methods are that their accuracy and stability depend on the time step size. The NSFD method guarantees a positive discrete solution for positive initial conditions. On the other hand, the disadvantage of the structured NSFD method is that for very large step sizes, a slight delay may occur in the traveling wave ^[30,33].

In this paper, NSFD method is used for discretization. Therefore, we need to provide the necessary information about this subject. NSFD method was introduced by Mickens's in 1989 ^[27]. Mickens shows that spurious solutions and numerical instabilities can be eliminated using this method. Mickens's use of the suitable denominator function played an important role in eliminating these problems and finding numerical solutions for each time-step size. For more information, see ^{[26,29,30,32]}.

The NSFD method for ordinary differential equations can be summarized as follows. We consider

where $\psi$ is a parameter. For this equation NSFD method is

In (3.3), $\Phi (h, \psi)$ is a denominator function and it can be chosen as: $\Phi (h, \psi) = \frac{1-e^{-\psi h}}{\psi }$. $\Phi (h, \psi)$ depends on the step size $h$ and $\psi$ is calculated from the information of fixed points of (3.2) ^[28]. The NSFD method can also be applied to fractional-order differential equations with the help of Grünwald-Letnikov discretization given in (2.8).

3.1.   Discretizations and stability analysis

In this discretization, we use the NSFD method, Grünwald-Letnikov definition for distributed order equations, and quadrature rule.

where $q_{0}^{\alpha_i} = (\Phi_k (h))^{-\alpha_i}, \ \ k = 1, 2, ..., 5$ and $0 < \alpha_i < 1$. The denominator functions are chosen as follows:

The expression $\sum_{j = 0}^{n+1} q_{j}^{\alpha_i} P_{n+1-j}$ in (3.4) can be expressed as:

If these adjustments are replaced in (3.4), $P_{n+1}$ is obtained by

If the same procedure is applied for (3.5)–(3.8) respectively, the discretized forms are obtained as follows:

This discrete form is expressed as system (3.14) for $\sum_{i = 1}^D \frac {w(\alpha_i)}{D} = K$ and $\sum_{i = 1}^D \frac {w(\alpha_i)}{D} (\Phi_j (h))^{-\alpha_i} = Z_j, j = 1, ..., 5$.

It is necessary to prove that all variables are not negative and $\tau_{2} < 1$ for system (3.14). Because in this case, the solutions of the system with positive initial conditions are positive for all $n > 0$.

Remark 3.1. If $P_{n} > 0, \ L_{n} > 0, \ S_{n} > 0, \ Q_{n} > 0, \ R_{n} > 0$ all variables are not negative ($\tau_{2} < 1$) and all the following conditions are satisfy, the solutions $P_{n+1}, \ L_{n+1}, \ S_{n+1}, \ Q_{n+1}$ and $R_{n+1}$ of system (3.14) are positive for all $n > 0$.

i) $\rho > K (q_{1}^{\alpha_i} P_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} P_{n+1-j}),$

ii) $\mu P_{n} S_{n} > K (q_{1}^{\alpha_i} L_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} L_{n+1-j}),$

iii) If $((Z_3)^{-\alpha_i}+\rho +\tau_{1} > \beta L_{n})$, $\gamma Q_{n} > K (q_{1}^{\alpha_i} S_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} S_{n+1-j})$ or

If $((Z_3)^{-\alpha_i}+\rho +\tau_{1} < \beta L_{n})$, $\gamma Q_{n} < K (q_{1}^{\alpha_i} S_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} S_{n+1-j}),$

iv) $\tau_{1} (1-\tau_{2}) S_{n} > K (q_{1}^{\alpha_i} Q_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} Q_{n+1-j}),$

v) $\tau_{2} \tau_{1} S_{n} > K (q_{1}^{\alpha_i} R_{n}+\sum_{j = 2}^{n+1} q_{j}^{\alpha_i} R_{n+1-j}).$

Now, the Jacobian $(J)$ matrix is required for the stability of this equilibrium point. So, the Jacobian matrix of system (3.14) is found as:

where,

Theorem 3.2. There are 2 types of equilibrium points for the smoking model (3.14), let's take $m_i = ((Z_i)^{-\alpha_i}+\rho+K \sum_{j = 1}^{n+1} q_{j}^{\alpha_i})$ for $i = 1, 2, ..., 5,$

i) If $S = 0$, smoking free equilibrium point $E = (\frac{\rho}{m_1}, 0, 0, 0, 0)$.

ii) If $S \ne 0$, system (3.14) has positive smoking present equilibrium point $E^* = (P^*, L^*, S^*, Q^*, R^*)$, where$P^* = \frac{ \rho}{(m_1+\mu S^*)}, \ L^* = \frac{ \mu (\frac{ \rho}{(m_1+\mu S^*)}) S^*}{(m_2+ \beta S^*)}, \ Q^* = \frac{ \tau_{1} (1-\tau_{2}) S^*}{(m_4 + \gamma)}, \ R^* = \frac{ \tau_{2} \tau_{1} S^* }{m_5},$ and $S^*$ is a positive solution of (3.16),

Proof. In order to find the equilibrium point of (3.14), the following statements should be provided.

Using (3.17) and $m_i = ((Z_i)^{-\alpha_i}+\rho+K \sum_{j = 1}^{n+1} q_{j}^{\alpha_i})$ for $i = 1, 2, ..., 5$, we obtain $P = \frac{\rho}{(m_1+\mu S)}$. Using the help of this expression and (3.18), we can find $L = \frac{ \mu (\frac{ \rho}{(m_1+\mu S)}) S}{(m_2+ \beta S)}$. Likewise, if $Q$ and $R$ expressions are left alone in (3.20) and (3.21), $Q = \frac{ \tau_{1} (1-\tau_{2}) S}{(m_4+ \gamma)}$ and $R = \frac{ \tau_{2} \tau_{1} S }{(m_5)}$ are found.

i) If we choose $S = 0$ in (3.17)–(3.21), we can find $E = (\frac{ \rho}{m_1}, 0, 0, 0, 0)$ equilibrium point.

ii) If $S \ne 0$, by substituting the found expressions in (3.19),

With this equation, the solution of $S^*$ can be obtained. Thus the equilibrium point for this situation can be expressed as:

Theorem 3.3. The smoking free equilibrium point $E = (\frac{ \rho}{m_1}, 0, 0, 0, 0)$ is locally asymptotically stable if $R_o < 1$ and the following conditions are satisfied, if not unstable.

i) $|- K (q_{1}^{\alpha_i})| < |((Z_1)^{-\alpha_i}+\rho)|$

ii) $|- K (q_{1}^{\alpha_i})| < |((Z_2)^{-\alpha_i}+\rho)|$

iii) $| - K (q_{1}^{\alpha_i})| < |((Z_5)^{-\alpha_i}+\rho)|$

iv) $|-a_1+\sqrt{a_1^2-4a_1 a_0}| < 2$

v) $|-a_1-\sqrt{a_1^2-4a_1 a_0}| < 2$

where

Proof. The Jacobian matrix is:

With the help of this matrix, characteristic equation calculated as:

The eigenvalues of $J(E)$ are $\lambda_1 = \frac{- K (q_{1}^{\alpha_i})}{((Z_1)^{-\alpha_i}+\rho)}$, $\lambda_2 = \frac{- K (q_{1}^{\alpha_i})}{((Z_2)^{-\alpha_i}+\rho)}$, $\lambda_3 = \frac{ - K (q_{1}^{\alpha_i})}{((Z_5)^{-\alpha_i}+\rho)}$, $\lambda_4 = \frac{-a_1+\sqrt{a_1^2-4a_1 a_0}}{2}$ and $\lambda_5 = \frac{-a_1-\sqrt{a_1^2-4a_1 a_0}}{2}$. We know that $E$ is stable for all eigenvalues are less than one in absolute values with Lemma (2.7). So if $(i) - (v)$ conditions are satisfied, the equilibrium point will be stable. And also, if $R_o < 1$ that means $\tau_1 (1-\tau_2) \gamma > (\rho+\tau_1)(\rho+\gamma)$, $a_1^2-4a_1 a_0$ is less than $0$. So, $\lambda_4$ and $\lambda_5$ expressions are obtained as complex roots.

Therefore, if $R_o < 1$, $E$ is locally asymptotically stable because all eigenvalues have negative real parts.

Remark 3.4. The stability of the equilibrium point $E^*$ involves quite complex operations. Therefore, the stability of the $E^*$ point has been investigated in the numerical analysis section with the Schur-Cohn test.

4.   Numerical simulations

In this section, we have given some numerical simulations. These simulations were obtained by the NSFD method. All initial conditions and parameters are listed in Table 2 ^[40]. It is seen that the conditions given by Remark (3.1) are satisfy with Table 2 values, so the results of (3.14) give positive solutions. These positive results are shown in Figures 1–6.

$\bf{i) \ Analysis \ \ of \ \ E \ \ equilibrium \ \ point:}$ With the help of the values given in Table 2, $w(\alpha) = \alpha-0.4$ and $h = 0.4$, firstly we calculate reproductive number

The eigenvalues are obtained as:

From Lemma (2.8), we get $a_1 = 0.0000078$ and $a_0 = 0.02$. We show that

i) $\ 1+a_0+a_1 = 1.0200078 > 0,$

ii) $\ 1+a_0-a_1 = 1.0199922 > 0,$

iii) $\ a_0 = 0.02 < 1,$

where $|\lambda_i| < 1, \ i = 1, 2, 3, 4, 5.$ So, $E$ equilibrium point is asymptotically stable as all the conditions in Lemma (2.8) are satisfied.

$\bf{ii) \ Analysis \ \ of \ \ E^* \ \ equilibrium \ \ point:}$ With the same procedure of analysis of $E^*$, reproductive number calculated is

Using Lemma (2.9),

i) $\ 1+a_0+a_1+a_2+a_3+a_4 = 10.8305 > 0,$

ii) $\ 1-a_0+a_1-a_2+a_3-a_4 = 7.0784 > 0,$

iii) $\ a_0 = 0.0002 < 1,$

iv) $\ |c_{2, 0}| = 1.0000 > |c_{2, 4}| = 0.1698,$

v) $\ |c_{3, 0}| = 1.0180 < |c_{3, 3}| = 6.0244,$

vi) $\ |c_{4, 0}| = 36.9762 > |c_{4, 2}| = 11.3052.$

As you can see, Lemma (2.9) doesn't satisfy because of the condition $(v)$. So according to the Schur-Cohn test, $E^*$ is unstable.

In Figures 1–5, graphics of all the definitions in the smoking model are given ($\alpha = 1, h = 0.01$). We can observe that the number of potential smokers increases over time and the number of occasional, heavy, temporary smokers, and permanent quitters decreases over time. In here, the value of $w(\alpha)$ was taken as $w1 = 1, \ w2 = \alpha+0.3, \ w3 = \alpha-0.4, \ w4 = 2\alpha-0.75$ and $w5 = \alpha+0.85$. We see the results are consistent when the graphs were compared by Singh et al.^[40] for $w(\alpha) = 1$. Thus, it was seen that we can interpret the selection of the density function in this way about the fractional-order state of the equation with the help of the distributed order equation. Although $w(\alpha)$ has different options, we see the solutions got by the NSFD method approach the correct endemic equilibrium point. We have found a general solution, including the fractional type solution obtained in Singh et al. ^[40]. For this reason, we can understand the dynamics of diseases that may occur in the face of different external factors with the solutions obtained.

Finally, we show some phase portraits of the equilibrium point $E$ in Figure 6. We can measure biologically the behavior of all subgroups with the selection of the density function when all graphs are considered. Thus, we can predict the behavior of populations belonging to the model under variable conditions. According to these predictions, we can take precautions against various diseases.

In Table 3, we present the effect of time step size on the Runge-Kutta method, Theta method, and NSFD method. Here, we see that the NSFD discretization is more effective than the classical method for bigger step-size. In Table 4, we compare CPU times for three numerical methods. We can say that the numerical methods evaluated among themselves have not been extremely different. Here, we carried out numerical calculations using MATLAB to illustrate the dynamics of the system. Because of this study and simulations, we have seen that the use of distributed order differential equations is quite suitable for this model. Because the selection of the density function provided a fast and accurate interpretation for different situations.

5.   Conclusions

In this paper, we have established a distributed order differential equation system for the smoking model. We have used the NSFD method and approximate Grünwald-Letnikov formula for the numerical solution of this model. Then, we have expressed the graphics of the solutions. We also analyzed the equilibrium points with the help of graphics and lemmas. We can see from the graphics that the interpretation capability is quite easy thanks to the different options of the density function. As can be understood from these data, the findings cover all the literature studies on this subject. Because of the density function, we can transform the distributed order differential equations into other equations. Since time is very important in such an important health problem, it is vital to use the solutions found, as they give the characteristics of the populations even in different situations. For these reasons, the use of distributed order differential equations is very convenient both in this model and in most models that cause such health problems.

Acknowledgements

The author Mehmet KOCABIYIK would like to thank the Scientific and Technological Research Council of Turkey (TUBITAK 2211-E program).

Conflict of interest

All authors declare no conflicts of interest in this paper.