# Numerical Solution of Fuzzy Fractional Pharmacokinetics Model Arising from Drug Assimilation into the Bloodstream

Authors: Ali Ahmadian; Norazak Senu; Farhad Larki; Soheil Salahshour; Mohamed Suleiman; Md. Shabiul Islam
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We propose a Jacobi tau method for solving a fuzzy fractional pharmacokinetics. This problem can model the concentration of the drug in the blood as time increases. The proposed approach is
based on the Jacobi tau (JT) method. To illustrate the reliability of the method, some special cases of the equations are solved as test examples. The method reduces the solution of the problem to the solution of
a system of algebraic equations. Error analysis included the fractional derivative error estimation, and the upper bound of the absolute errors is introduced for this method.


---

## Body

## 1. Introduction

Pharmacokinetics is defined as the study of the time course of drug absorption, distribution, metabolism, and excretion [1]. Pharmacokinetics determines bioavailability, volume of distribution, and clearance. Bioavailability is the fraction of a drug absorbed into the systemic circulation, while volume of distribution and clearance are measure of apparent space in the body available to contain the drug and measure of the body’s ability to eliminate the drug, respectively. In pharmacokinetics, as a basic hypothesis, always a relationship exists between the pharmacologic or toxic response to a drug and the concentration of the drug in the blood (plasma). However, for some drugs there is no straightforward relationship between concentration in plasma and pharmacologic effects [2]. Generally, the concentration of drug in the systemic circulation is related to the concentration of drug at its sites of action. Absorption, drug distribution (target tissues), drug concentration in circulation system, rate of elimination/plasma concentration, and elimination (metabolism, excretion) are some parameters that directly modify the pharmacokinetics parameters.Mathematics is widely used for the quantitative description of drug absorption, distribution, metabolism, and excretion (ADME). Some parameters in pharmacokinetics are defined by mathematic equations which can be obtained directly by measurement or through calculation using experimental data based on developed mathematical equations. For example, the first and simplest model used for drug absorption from the gastrointestinal (GI) tract assumes a “pseudosteady state” and utilizes the physicochemical properties of the compound in conjunction with the “pH-partition” hypothesis to predict the fraction of dose absorbed [3–5]. Tissue distribution is another important determinant of the pharmacokinetics profile of a drug. Hence, in drug development, the prediction of tissue distribution would help predict the in vivo pharmacokinetics of a compound prior to any experiments in animals or man. Two methods are available to predict tissue distribution. These methods predict either tissue: plasma ratios or the volume of distribution at steady state (V). From molecular descriptors, properties such as lipophilicity/hydrophobicity are estimated [6]. In some approaches, plasma protein binding is also taken into consideration [7]. Physiological information on tissue composition, the blood composition, and blood flow to the tissue is utilized to develop a partitioning model [8–10]. Möhler et al. [11] proposed a model that describes cell growth and detailed investigations on the metabolism and the kinetics of the influenza infection cycle which allow for the optimization of influenza virus vaccine production. In addition to these simple models, more complex mathematical manipulations called mathematical models have also been used to describe pharmacokinetics [12]. There are many classifications of pharmacokinetics models. In a specific classification for pharmacokinetics modeling, two general approaches can be considered: compartment-based modeling [13] and noncompartment-based modeling [14]. Both types of modeling take advantage of the quantitative structure-pharmacokinetics relationships that are described by empirical mathematical algorithms. They can be used to estimate the activity of a compound based on its chemical structure in a numeric format [15].In this study, our focus is to find the approximate solution of the fuzzy fractional model of the compartment model for the flow of antihistamine in the blood. To this end, we use the operational matrix of the Caputo fractional derivative of the Jacobi tau approximation based on the Jacobi polynomials to derive the fuzzy fractional approximate solution under Hukuhara differentiability (H-differentiability). To the best of our knowledge, there are no results on the JT method for solving pharmacokinetics equations arising in mathematical physics. This partially motivated our interest in such model and solution method.


### 1.1. Compartment Model

The compartment model framework is an extremely natural and valuable means with which to formulate models for processes which have inputs and outputs over time [16]. Compartment model itself can be classified in one-compartment model and two-compartment model. In the former one, the body is depicted as a kinetically homogeneous unit, while in the latter one the body resolves into a central compartment and a peripheral compartment [17]. According to the one-compartment model, the drug instantaneously distributes throughout the body, and at the same time it equilibrates between tissues. Thus, the drug concentration time profile shows a monophasic response. On the other hand, in two-compartment model the drug does not achieve instantaneous equilibration between the two compartments. Although, in this model, two compartments have no anatomical or physiological meaning, it is usually surmised that the central compartment is composed of tissues that are highly perfused such as heart, lungs, kidneys, liver, and brain while the peripheral compartment comprises less well-perfused tissues such as muscle, fat, and skin [17].


### 1.2. Methods Have Been Proposed to Solve Fractional Pharmacokinetics Models

While an increasing number of fractional order integrals and differential equations applications have been reported in the physics [18, 19], signal processing [20], engineering, and bioengineering literatures [21, 22], little attention has been paid to this class of models in the pharmacokinetics-pharmacodynamic (PKPD) literature. The reasons for the lack of application to pharmacodynamics are mainly two: first and foremost, PKPD models incorporating fractional calculus have not been proposed; second is computational: while the analytical solution of fractional differential equations is available in special cases, it turns out that even the simplest PKPD models that can be constructed using fractional calculus do not allow analytical solutions [23, 24].In this paper, we propose new category of PKPD models incorporating fuzzy fractional calculus and investigate their behavior using purposely written algorithm. The main purpose of this paper is to attract the attention of the field into these possibly interesting families of PKPD models. In particular, we will focus on compartment models, direct concentration-response relationships describing models relating drug concentration (in plasma or biophase) to a pharmacodynamic effect.The concept of fractional or noninteger order derivation and integration can be traced back to the genesis of integer order calculus itself [25–27]. Due to its tremendous scope and applications in several disciplines, a considerable attention has been given to exact and numerical solutions of fractional differential equations which is an extremely difficult task. Moreover, the solution techniques and their reliability are rather more important aspects. Several methods have been proposed by many researchers to solve the fractional order differential equations such as Taylor series [28], variational iteration method [29–31], adomian decomposition method [32, 33], fractional differential transform method (FDTM) [34, 35], and homotopy analysis method [36, 37].As it is known, the spectral method is one of the flexible methods of discretization for most types of differential equations [38–40]. Historically, spectral method has been relegated to fractional calculus, but in few years, it has been successfully applied for the fractional equation models based on the different types of orthogonal polynomials such as Block pulse functions [41, 42], Legendre polynomials [43–46], Chebyshev polynomials [47–49], Laguerre polynomials [50–53], and Bernstein polynomials [54–56]. Doha et al. [57] introduced the shifted Jacobi operational matrix of fractional derivative which is based on Jacobi tau method for solving numerically linear multiterm fractional differential equations with initial or boundary conditions. Subsequently, Kazem [58] generalized Jacobi integral operational matrix to fractional calculus. Afterwards, Doha et al. [59] proposed a direct solution techniques for solving the linear multiorder fractional initial value problems with constant and variable coefficients using shifted Jacobi tau method and quadrature shifted Jacobi tau method, respectively. In this paper, we intend to extend the application of the Jacobi polynomials to solve fuzzy fractional pharmacokinetics model of order[0,1].On the other hand, a considerable attention has been made to the fractional differential equations in the sense of fuzzy setting theory. Agarwal et al. [60] is among the pioneers who presented the concept of solutions for fractional differential equations with uncertainty. Thereafter, the existence and uniqueness of the solution of the fuzzy fractional differential equations (FFDEs) under the Riemann-Liouville and Caputo’s fuzzy fractional differentiability were investigated in the literature [61–64]. However, the application of the numerical methods for solving FFDEs is unknown and traceless in the literature of the fuzzy fractional calculus, and only a few number of researches have been reported for the approximate solution of FFDEs [64–69].In the present paper, we intend to introduce new families of PKPD models based on the application of fractional calculus to PKPD models. The aim of the paper is not to claim the superiority of fractional dynamics models with respect to standard ones, but it is simply to define the new families and provide some insights into their qualitative behavior. The main purpose is to apply the fuzzy logic in differential equations of fractional order which has been used as an effective tool for considering uncertainty in modeling the processes; FFDEs can also offer a more comprehensive account of the process or phenomenon, specifically for analyzing the behavior of the PKPD models. Furthermore, we suggest possible applications and stimulate further research, which might, or might not, demonstrate the applicability and importance of spectral methods by using of the orthogonal polynomials for finding the approximate solutions of the PKPD models based on the fuzzy fractional calculus.The structure of this paper is as follows. In the next section, we briefly recall the mathematical foundations of fractional calculus, required definitions of fuzzy setting theory, and summarize the properties of Jacobi polynomials. In Section3, we then provide illustrations of the governing fraction equation. In Section 4, the proposed method is explained for numerical solution of the derived FFDE. Section 5 is devoted to the numerical solution of the problem with different dose of the drug, and the error analysis is provided to demonstrate the applicability and validity of the method. A final discussion comments on the results are presented in Section 6.


## 1.1. Compartment Model

The compartment model framework is an extremely natural and valuable means with which to formulate models for processes which have inputs and outputs over time [16]. Compartment model itself can be classified in one-compartment model and two-compartment model. In the former one, the body is depicted as a kinetically homogeneous unit, while in the latter one the body resolves into a central compartment and a peripheral compartment [17]. According to the one-compartment model, the drug instantaneously distributes throughout the body, and at the same time it equilibrates between tissues. Thus, the drug concentration time profile shows a monophasic response. On the other hand, in two-compartment model the drug does not achieve instantaneous equilibration between the two compartments. Although, in this model, two compartments have no anatomical or physiological meaning, it is usually surmised that the central compartment is composed of tissues that are highly perfused such as heart, lungs, kidneys, liver, and brain while the peripheral compartment comprises less well-perfused tissues such as muscle, fat, and skin [17].


## 1.2. Methods Have Been Proposed to Solve Fractional Pharmacokinetics Models

While an increasing number of fractional order integrals and differential equations applications have been reported in the physics [18, 19], signal processing [20], engineering, and bioengineering literatures [21, 22], little attention has been paid to this class of models in the pharmacokinetics-pharmacodynamic (PKPD) literature. The reasons for the lack of application to pharmacodynamics are mainly two: first and foremost, PKPD models incorporating fractional calculus have not been proposed; second is computational: while the analytical solution of fractional differential equations is available in special cases, it turns out that even the simplest PKPD models that can be constructed using fractional calculus do not allow analytical solutions [23, 24].In this paper, we propose new category of PKPD models incorporating fuzzy fractional calculus and investigate their behavior using purposely written algorithm. The main purpose of this paper is to attract the attention of the field into these possibly interesting families of PKPD models. In particular, we will focus on compartment models, direct concentration-response relationships describing models relating drug concentration (in plasma or biophase) to a pharmacodynamic effect.The concept of fractional or noninteger order derivation and integration can be traced back to the genesis of integer order calculus itself [25–27]. Due to its tremendous scope and applications in several disciplines, a considerable attention has been given to exact and numerical solutions of fractional differential equations which is an extremely difficult task. Moreover, the solution techniques and their reliability are rather more important aspects. Several methods have been proposed by many researchers to solve the fractional order differential equations such as Taylor series [28], variational iteration method [29–31], adomian decomposition method [32, 33], fractional differential transform method (FDTM) [34, 35], and homotopy analysis method [36, 37].As it is known, the spectral method is one of the flexible methods of discretization for most types of differential equations [38–40]. Historically, spectral method has been relegated to fractional calculus, but in few years, it has been successfully applied for the fractional equation models based on the different types of orthogonal polynomials such as Block pulse functions [41, 42], Legendre polynomials [43–46], Chebyshev polynomials [47–49], Laguerre polynomials [50–53], and Bernstein polynomials [54–56]. Doha et al. [57] introduced the shifted Jacobi operational matrix of fractional derivative which is based on Jacobi tau method for solving numerically linear multiterm fractional differential equations with initial or boundary conditions. Subsequently, Kazem [58] generalized Jacobi integral operational matrix to fractional calculus. Afterwards, Doha et al. [59] proposed a direct solution techniques for solving the linear multiorder fractional initial value problems with constant and variable coefficients using shifted Jacobi tau method and quadrature shifted Jacobi tau method, respectively. In this paper, we intend to extend the application of the Jacobi polynomials to solve fuzzy fractional pharmacokinetics model of order[0,1].On the other hand, a considerable attention has been made to the fractional differential equations in the sense of fuzzy setting theory. Agarwal et al. [60] is among the pioneers who presented the concept of solutions for fractional differential equations with uncertainty. Thereafter, the existence and uniqueness of the solution of the fuzzy fractional differential equations (FFDEs) under the Riemann-Liouville and Caputo’s fuzzy fractional differentiability were investigated in the literature [61–64]. However, the application of the numerical methods for solving FFDEs is unknown and traceless in the literature of the fuzzy fractional calculus, and only a few number of researches have been reported for the approximate solution of FFDEs [64–69].In the present paper, we intend to introduce new families of PKPD models based on the application of fractional calculus to PKPD models. The aim of the paper is not to claim the superiority of fractional dynamics models with respect to standard ones, but it is simply to define the new families and provide some insights into their qualitative behavior. The main purpose is to apply the fuzzy logic in differential equations of fractional order which has been used as an effective tool for considering uncertainty in modeling the processes; FFDEs can also offer a more comprehensive account of the process or phenomenon, specifically for analyzing the behavior of the PKPD models. Furthermore, we suggest possible applications and stimulate further research, which might, or might not, demonstrate the applicability and importance of spectral methods by using of the orthogonal polynomials for finding the approximate solutions of the PKPD models based on the fuzzy fractional calculus.The structure of this paper is as follows. In the next section, we briefly recall the mathematical foundations of fractional calculus, required definitions of fuzzy setting theory, and summarize the properties of Jacobi polynomials. In Section3, we then provide illustrations of the governing fraction equation. In Section 4, the proposed method is explained for numerical solution of the derived FFDE. Section 5 is devoted to the numerical solution of the problem with different dose of the drug, and the error analysis is provided to demonstrate the applicability and validity of the method. A final discussion comments on the results are presented in Section 6.


## 2. Preliminaries and Notations

In this section, we are going to state the definition and preliminaries of fuzzy mathematics [63, 70, 71], fractional calculus [25, 27], and some properties of shifted Jacobi polynomials.


### 2.1. The Fuzzy Settings Definitions

Definition 1.
Letube a fuzzy set inℝ.uis called a fuzzy number if:(i)
uis normal: there existsx0∈ℝsuch thatu(x0)=1;(ii)
uis convex: for allx,y∈ℝand0≤λ≤1, it holds that
(1)u(λx+(1-λ)y≥min{u(x),u(y)});(iii)
uis upper semicontinuous: for anyx0∈ℝ, it holds that
(2)u(x0)≥limx→x0±u(x);(iv)
[u]0=supp(u)¯is a compact subset ofℝ.
In this paper, the set of all fuzzy numbers is denoted byℝℱ.Definition 2.
Letu∈ℝℱandr∈[0,1]. Ther-cut ofuis the crisp set[u]rthat contains all elements with membership degree inugreater than or equal tor; that is
(3)[u]r={x∈ℝ∣u(x)≥r}.
For a fuzzy numberu, itsr-cuts are closed intervals inℝ, and we denote them by
(4)[u]r=[u1r,u2r].According to Zadeh’s extension principle, the operation of addition onℝℱis defined as follows:
(5)(u+v)(x)=supy∈ℝmin{u(y),v(x-y)},x∈ℝ,
and scalar multiplication of a fuzzy number is given by
(6)(k⊙u)(x)={u(xk),k>0,0~,k=0,
where0~∈ℝℱ.Definition 3 (see [70]).
The distanceD(u,v)between two fuzzy numbersuandvis defined as
(7)D(u,v)=supr∈[0,1]dH([u]r,[v]r),
where
(8)dH([u]r,[v]r)=max{|u1r-v1r|,|u2r-v2r|}
is the Hausdorff distance between[u]rand[v]r.It is easy to see thatDis a metric inℝℱand has the following properties (see [70, 72]):(i)
D(u⊕w,v⊕w)=D(u,v),∀u,v,w∈ℝℱ,(ii)
D(k⊙u,k⊙v)=|k|D(u,v),∀k∈ℝ,u,v∈ℝℱ,(iii)
D(u⊕v,w⊕e)≤D(u,w)+D(v,e),∀u,v,w∈ℝℱ,(iv)
D(u+v,0)≤D(u,0)+D(v,0),∀u,v∈ℝℱ,(v)
(ℝℱ,D)is a complete metric space.Definition 4 (see [73]).
Letfandgbe the two fuzzy-number-valued functions on the interval[a,b],that is,f,g:[a,b]→ℝℱ. The uniform distance between fuzzy-number-valued functions is defined by
(9)D*(f,g):=supx∈[a,b]D(f(x),g(x)).Remark 5 (see [73]).
Letf:[a,b]→ℝℱbe fuzzy continuous. Then from property (iv) of Hausdorff distance, we can define
(10)D(f(x),0~)=supr∈[0,1]max{|f1r(x)|,|f2r(x)|},∀x∈[a,b].Definition 6 (see [74]).
Letx,y∈ℝℱ. If there existsz∈ℝℱsuch thatx=y⊕z, thenzis called the H-difference ofxandy, and it is denoted byx⊖y.In this paper, the sign “⊖” always stands for H-difference, and note thatx⊕y≠x+(-y). Also throughout the paper, it is assumed that the Hukuhara difference and generalized Hukuhara differentiability exist.Theorem 7 (see [75]).
LetF:(a,b)→ℝℱbe a function and denote[F(t)]r=[fr(t),gr(t)], for eachr∈[0,1]. Then,(1)
ifFis (1)-differentiable, thenfr(t)andgr(t)are differentiable functions and
(11)[F′(t)]r=[fr′(t),gr′(t)];(2)
ifFis (2)-differentiable, thenfr(t)andgr(t)are differentiable functions and
(12)[F′(t)]r=[gr′(t),fr′(t)].Definition 8 (see [76]).
Consider then×nlinear system of the following equations:
(13)a11x1+a12x2+⋯+a1nxn=y1,a21x1+a22x2+⋯+a2nxn=y2,⋮an1x1+an2x2+⋯+annxn=yn.
The matrix form of the above equations is
(14)AX=Y,
where the coefficient matrixA=(aij),1≤i,j≤nis a crispn×nmatrix andyi∈ℝℱ,1≤i≤n. This system is called a fuzzy linear system (FLS).Definition 9 (see [76]).
A fuzzy number vector(x1,x2,…,xn)tgiven byxi=(xi-r,xi+r),1≤i≤n,0≤r≤1is called a solution of the fuzzy linear system (2) if
(15)(∑j=1n‍aijxj)-r=∑j=1n‍(aijxj)-r=yi-r,(∑j=1n‍aijxj)+r=∑j=1n‍(aijxj)+r=yi+r.
If for a particulark,akj>0,1≤j≤n, we simply get
(16)∑j=1n‍akjxj-r=yk-r,∑j=1n‍akjxj+r=yk+r.
To solve fuzzy linear systems, see [77].In this part, we firstly give some basic definitions and some properties of fractional calculus [27]. Afterwards, the extension of the fractional differentiability in the sense of fuzzy concept is provided, and some relevant properties which are used in the rest of the paper are given [61, 63, 64].Letmbe the smallest integer that exceedsv, then Caputo’s fractional derivative operator of order v>0is defined as
(17)cDvf(x)={Jm-vDmf(x),ifm-1<v<m,Dmf(x),ifv=m,m∈ℕ,
where
(18)Jvf(x)=1Γ(v)∫0x(x-t)v-1f(t)dt,v>0,x>0.
For the Caputo derivative, we have
(19)cDvC=0,(Cisaconstant),cDvxβ={0,forβ∈ℕ0,β<⌈v⌉,Γ(β+1)Γ(β+1-v)xβ-v,forβ∈ℕ0,β≥⌈v⌉orβ∉ℕ,β>⌊v⌋.
Caputo’s fractional differentiation is a linear operation, namely,
(20)cDv(λf(x)+μg(x))=λcDvf(x)+μcDvg(x),
whereλandμare constants.Theorem 10 (generalized Taylor formula [78]).
Suppose thatDakαf(t)∈C(a,b]fork=0,1,…,n+1where0<α≤1; then
(21)f(t)=∑i=0n‍(t-a)iαΓ(iα+1)[cDakαf(t)]t=a+Rnα(t,a),
with
(22)Rnα(t,a)=(t-a)(n+1)αΓ((n+1)α+1)[cDa(n+1)αf(t)]t=ξ,a⩽ξ⩽t,t∈[a,b],
where
(23)cDanα=cDaαccDaα,…,cDaα︸ntimes.Here,cDαshows fractional derivative operator in the Caputo sense, and fTni(t)=∑i=0n((t-a)iα/Γ(iα+1))[cDakαf(t)]t=ais called fractional Taylor expansion with degreenioff, andRnα(t,a)is reminder term in fractional Taylor expansion off.Now, we present some fuzzy fractional notations which are used later in the paper.(i)
Lpℝℱ(a,b),1≤p≤∞is the set of all fuzzy-valued measurable functionsfon[a,b]where||f||p=(∫01‍(d(f(t),0))pdt)1/p.(ii)
Cℝℱ[a,b]is a space of fuzzy-valued functions which are continuous on[a,b].(iii)
Cnℝℱ[a,b]indicates the set of all fuzzy-valued functions which are continuous up to ordern.(iv)
ACℝℱ[a,b]denotes the set of all fuzzy-valued functions which are absolutely continuous.Note that one can easily find these notations in the crisp context in [25, 27] and references therein.Definition 11 (see [63]).
Letf∈Cℝℱ[a,b]∩Lℝℱ[a,b]. The Riemann-Liouville integral of fuzzy-valued functionfis defined as
(24)(RLIa+vf)(x)=1Γ(v)∫ax‍f(t)dt(x-t)1-v,x>a,0<v≤1.Definition 12 (see [63]).
Letf∈Cℝℱ[a,b]∩Lℝℱ[a,b]. Then, f is said to be Caputo’s H-differentiable atxwhen
(25)(i)(cDa+vf)(x)=(RLDa+β[f(t)⊖f(a)])(x),(ii)(cDa+vf)(x)=(RLDa+β[-f(a)⊖(-f(t))])(x).Definition 13 (see [63]).
Letf:Lℝℱ[a,b]∩Cℝℱ[a,b]andx0∈(a,b)and Φ(x)=(1/Γ(1-v))∫ax‍(f(t)/(x-t)v)dt. We say that f(x)is fuzzy Caputo fractional differentiable of order 0<v≤1atx0, if there exists an element(cDa+vf)(x0)∈Cℝℱ[a,b][a,b]such that for all0≤r≤1,h>0,
(26)(i)(cDa+vf)(x0)=limh→0+Φ(x0+h)⊖Φ(x0)h=limh→0+Φ(x0)⊖Φ(x0-h)h,
or
(27)(ii)(cDa+vf)(x0)=limh→0+Φ(x0)⊖Φ(x0+h)-h=limh→0+Φ(x0-h)⊖Φ(x0)-h,
or
(28)(iii)(cDa+vf)(x0)=limh→0+Φ(x0+h)⊖Φ(x0)h=limh→0+Φ(x0-h)⊖Φ(x0)-h,
or
(29)(iv)(cDa+vf)(x0)=limh→0+Φ(x0)⊖Φ(x0+h)-h=limh→0+Φ(x0)⊖Φ(x0-h)h.For the sake of simplicity, we say that the fuzzy-valued functionfisc[(1)-v]-differentiable if it is differentiable as in Definition 13 case (i), andfisc[(2)-v]-differentiable if it is differentiable as in Definition 13 case (ii), and so on for the other cases.Theorem 14 (see [63]).
Let0<v≤1andf∈ACℝℱ[a,b]; then the fuzzy Caputo fractional derivative exists almost everywhere on(a,b), and for all0≤r≤1one has
(30)(cDa+vf)(x;r)=[1Γ(1-v)∫axf-′r(t)dt(x-t)v,1Γ(1-v)∫ax‍f+′r(t)dt(x-t)v]=[(Ia+1-vDf-r)(x),(Ia+1-vDf+r)(x)],
whenfis (1)-differentiable, and
(31)(cDa+vf)(x;r)=[1Γ(1-v)∫ax‍f+′r(t)dt(x-t)v,1Γ(1-v)∫ax‍f-′r(t)dt(x-t)v,]=[(Ia+1-vDf+r)(x),(Ia+1-vDf-r)(x)],
when f is (2)-differentiable.Theorem 15 (fuzzy generalized Taylor’s formula [72]).
Letf(x)∈ACℝℱ[a,b](0,b], and suppose that cDkαf(x)∈Cℝℱ[a,b](0,b]fork=0,1,…,n+1where0<α<1,0≤x0≤xandx∈(0,b]. Then, one has
(32)[f(x)]r=[f_r(x),f¯r(x)],f_r(x)=∑i=0nxiαΓ(iα+1)cDiαf_r(0+)+cD(n+1)αfr(x0)Γ(nα+α+1)x(n+1)α,f¯r(x)=∑i=0nxiαΓ(iα+1)cDiαf¯r(0+)+cD(n+1)αfr(x0)Γ(nα+α+1)x(n+1)α,
where cDαf_r(0)=cDαf_r(x)|x=0,cDαf¯r(0)=cDαf¯r(x)|x=0.


### 2.2. Jacobi Polynomials

The well-known Jacobi polynomials associated with the parameters (α>-1,β>-1) (see, e.g., Luke [79] and Szegö [80]), are a sequence of polynomialsPi(α,β)(t)(i=0,1,…), each, respectively, of degreei. For using these polynomials on(0,L), we present the shifted Jacobi polynomials by implementing the change of variablet=(2x/L-1). Let the shifted Jacobi polynomials Pi(α,β)(2x/L-1)be denoted byPL,i(α,β)(x), satisfying the orthogonality relation
(33)∫0L‍PL,j(α,β)(x)PL,k(α,β)(x)wL(α,β)(x)dx=hk,
wherewL(α,β)(x)=xβ(L-x)αand
(34)hk={Lα+β+1Γ(k+α+1)Γ(k+β+1)(2k+α+β+1)k!Γ(k+α+β+1)i=j,0,i≠j.
The shifted Jacobi polynomialPL,i(α,β)(x)of degreeihas the form
(35)PL,i(α,β)(x)=∑k=0i‍(-1)i-kΓ(i+β+1)Γ(i+k+α+β+1)Γ(k+β+1)Γ(i+α+β+1)(i-k)!k!Lkxk,
where
(36)PL,i(α,β)(0)=(-1)iΓ(i+β+1)Γ(β+1)i!,PL,i(α,β)(L)=Γ(i+α+1)Γ(α+1)i!.
Also we can state the shifted Jacobi polynomial by the following recurrence relation:
(37)PL,i(α,β)(x)=(α+β+2i-1){(α2-β2+(2xL-1))×(α+β+2i)(α+β+2i-2)(2xL-1)}×(2i(α+β+i)(α+β+2i-2))-1×PL,i-1(α,β)(x)-(α+i-1)(β+i-1)(α+β+2i)i(α+β+i)(α+β+2i-2)×PL,i-2(α,β)(x)i=2,3,…,
wherePL,0(α,β)(x)=1andPL,1(α,β)(x)=((α+β+2)/2)(2x/L-1)+(α-β)/2.We notice that a functionu(x), square integrable in(0,L), can be expanded in terms of shifted Jacobi polynomials as
(38)u(x)=∑j=0+∞‍ajPL,j(α,β)(x),
where the coefficientsajare
(39)aj=1hj∫0L‍PL,j(α,β)(x)u(x)wL(α,β)(x)dx,j=0,1,…
For a given particular problem, only the first(N+1)-terms shifted Jacobi polynomials are used. Therefore, we have
(40)uN(x)≃∑j=0N‍ajPL,j(α,β)(x).Theorem 16 (see [59]).
The Caputo fractional derivative of ordervof the shifted Jacobi polynomials of degreeiis obtained from
(41)DvPL,i(α,β)(x)=∑j=0∞‍Sv(i,j,α,β)PL,i(α,β)(x),i=⌈v⌉,⌈v⌉+1,…,
where
(42)Sv(i,j,α,β)=∑k=⌈v⌉i(-1)i-kLα+β-v+1Γ(j+β+1)×Γ(i+β+1)Γ(i+k+α+β+1)×(hjΓ(j+k+α+β+1)Γ(k+β+1)×Γ(i+α+β+1)Γ(k-v+1)(i-k)!)-1×∑l=0j(-1)j-lΓ(j+l+α+β+1)Γ(α+1)×Γ(l+k+β-v+1)×(Γ(l+β+1)Γ(l+k+α+β-v+2)(j-l)!l!)-1.


## 2.1. The Fuzzy Settings Definitions

Definition 1.
Letube a fuzzy set inℝ.uis called a fuzzy number if:(i)
uis normal: there existsx0∈ℝsuch thatu(x0)=1;(ii)
uis convex: for allx,y∈ℝand0≤λ≤1, it holds that
(1)u(λx+(1-λ)y≥min{u(x),u(y)});(iii)
uis upper semicontinuous: for anyx0∈ℝ, it holds that
(2)u(x0)≥limx→x0±u(x);(iv)
[u]0=supp(u)¯is a compact subset ofℝ.
In this paper, the set of all fuzzy numbers is denoted byℝℱ.Definition 2.
Letu∈ℝℱandr∈[0,1]. Ther-cut ofuis the crisp set[u]rthat contains all elements with membership degree inugreater than or equal tor; that is
(3)[u]r={x∈ℝ∣u(x)≥r}.
For a fuzzy numberu, itsr-cuts are closed intervals inℝ, and we denote them by
(4)[u]r=[u1r,u2r].According to Zadeh’s extension principle, the operation of addition onℝℱis defined as follows:
(5)(u+v)(x)=supy∈ℝmin{u(y),v(x-y)},x∈ℝ,
and scalar multiplication of a fuzzy number is given by
(6)(k⊙u)(x)={u(xk),k>0,0~,k=0,
where0~∈ℝℱ.Definition 3 (see [70]).
The distanceD(u,v)between two fuzzy numbersuandvis defined as
(7)D(u,v)=supr∈[0,1]dH([u]r,[v]r),
where
(8)dH([u]r,[v]r)=max{|u1r-v1r|,|u2r-v2r|}
is the Hausdorff distance between[u]rand[v]r.It is easy to see thatDis a metric inℝℱand has the following properties (see [70, 72]):(i)
D(u⊕w,v⊕w)=D(u,v),∀u,v,w∈ℝℱ,(ii)
D(k⊙u,k⊙v)=|k|D(u,v),∀k∈ℝ,u,v∈ℝℱ,(iii)
D(u⊕v,w⊕e)≤D(u,w)+D(v,e),∀u,v,w∈ℝℱ,(iv)
D(u+v,0)≤D(u,0)+D(v,0),∀u,v∈ℝℱ,(v)
(ℝℱ,D)is a complete metric space.Definition 4 (see [73]).
Letfandgbe the two fuzzy-number-valued functions on the interval[a,b],that is,f,g:[a,b]→ℝℱ. The uniform distance between fuzzy-number-valued functions is defined by
(9)D*(f,g):=supx∈[a,b]D(f(x),g(x)).Remark 5 (see [73]).
Letf:[a,b]→ℝℱbe fuzzy continuous. Then from property (iv) of Hausdorff distance, we can define
(10)D(f(x),0~)=supr∈[0,1]max{|f1r(x)|,|f2r(x)|},∀x∈[a,b].Definition 6 (see [74]).
Letx,y∈ℝℱ. If there existsz∈ℝℱsuch thatx=y⊕z, thenzis called the H-difference ofxandy, and it is denoted byx⊖y.In this paper, the sign “⊖” always stands for H-difference, and note thatx⊕y≠x+(-y). Also throughout the paper, it is assumed that the Hukuhara difference and generalized Hukuhara differentiability exist.Theorem 7 (see [75]).
LetF:(a,b)→ℝℱbe a function and denote[F(t)]r=[fr(t),gr(t)], for eachr∈[0,1]. Then,(1)
ifFis (1)-differentiable, thenfr(t)andgr(t)are differentiable functions and
(11)[F′(t)]r=[fr′(t),gr′(t)];(2)
ifFis (2)-differentiable, thenfr(t)andgr(t)are differentiable functions and
(12)[F′(t)]r=[gr′(t),fr′(t)].Definition 8 (see [76]).
Consider then×nlinear system of the following equations:
(13)a11x1+a12x2+⋯+a1nxn=y1,a21x1+a22x2+⋯+a2nxn=y2,⋮an1x1+an2x2+⋯+annxn=yn.
The matrix form of the above equations is
(14)AX=Y,
where the coefficient matrixA=(aij),1≤i,j≤nis a crispn×nmatrix andyi∈ℝℱ,1≤i≤n. This system is called a fuzzy linear system (FLS).Definition 9 (see [76]).
A fuzzy number vector(x1,x2,…,xn)tgiven byxi=(xi-r,xi+r),1≤i≤n,0≤r≤1is called a solution of the fuzzy linear system (2) if
(15)(∑j=1n‍aijxj)-r=∑j=1n‍(aijxj)-r=yi-r,(∑j=1n‍aijxj)+r=∑j=1n‍(aijxj)+r=yi+r.
If for a particulark,akj>0,1≤j≤n, we simply get
(16)∑j=1n‍akjxj-r=yk-r,∑j=1n‍akjxj+r=yk+r.
To solve fuzzy linear systems, see [77].In this part, we firstly give some basic definitions and some properties of fractional calculus [27]. Afterwards, the extension of the fractional differentiability in the sense of fuzzy concept is provided, and some relevant properties which are used in the rest of the paper are given [61, 63, 64].Letmbe the smallest integer that exceedsv, then Caputo’s fractional derivative operator of order v>0is defined as
(17)cDvf(x)={Jm-vDmf(x),ifm-1<v<m,Dmf(x),ifv=m,m∈ℕ,
where
(18)Jvf(x)=1Γ(v)∫0x(x-t)v-1f(t)dt,v>0,x>0.
For the Caputo derivative, we have
(19)cDvC=0,(Cisaconstant),cDvxβ={0,forβ∈ℕ0,β<⌈v⌉,Γ(β+1)Γ(β+1-v)xβ-v,forβ∈ℕ0,β≥⌈v⌉orβ∉ℕ,β>⌊v⌋.
Caputo’s fractional differentiation is a linear operation, namely,
(20)cDv(λf(x)+μg(x))=λcDvf(x)+μcDvg(x),
whereλandμare constants.Theorem 10 (generalized Taylor formula [78]).
Suppose thatDakαf(t)∈C(a,b]fork=0,1,…,n+1where0<α≤1; then
(21)f(t)=∑i=0n‍(t-a)iαΓ(iα+1)[cDakαf(t)]t=a+Rnα(t,a),
with
(22)Rnα(t,a)=(t-a)(n+1)αΓ((n+1)α+1)[cDa(n+1)αf(t)]t=ξ,a⩽ξ⩽t,t∈[a,b],
where
(23)cDanα=cDaαccDaα,…,cDaα︸ntimes.Here,cDαshows fractional derivative operator in the Caputo sense, and fTni(t)=∑i=0n((t-a)iα/Γ(iα+1))[cDakαf(t)]t=ais called fractional Taylor expansion with degreenioff, andRnα(t,a)is reminder term in fractional Taylor expansion off.Now, we present some fuzzy fractional notations which are used later in the paper.(i)
Lpℝℱ(a,b),1≤p≤∞is the set of all fuzzy-valued measurable functionsfon[a,b]where||f||p=(∫01‍(d(f(t),0))pdt)1/p.(ii)
Cℝℱ[a,b]is a space of fuzzy-valued functions which are continuous on[a,b].(iii)
Cnℝℱ[a,b]indicates the set of all fuzzy-valued functions which are continuous up to ordern.(iv)
ACℝℱ[a,b]denotes the set of all fuzzy-valued functions which are absolutely continuous.Note that one can easily find these notations in the crisp context in [25, 27] and references therein.Definition 11 (see [63]).
Letf∈Cℝℱ[a,b]∩Lℝℱ[a,b]. The Riemann-Liouville integral of fuzzy-valued functionfis defined as
(24)(RLIa+vf)(x)=1Γ(v)∫ax‍f(t)dt(x-t)1-v,x>a,0<v≤1.Definition 12 (see [63]).
Letf∈Cℝℱ[a,b]∩Lℝℱ[a,b]. Then, f is said to be Caputo’s H-differentiable atxwhen
(25)(i)(cDa+vf)(x)=(RLDa+β[f(t)⊖f(a)])(x),(ii)(cDa+vf)(x)=(RLDa+β[-f(a)⊖(-f(t))])(x).Definition 13 (see [63]).
Letf:Lℝℱ[a,b]∩Cℝℱ[a,b]andx0∈(a,b)and Φ(x)=(1/Γ(1-v))∫ax‍(f(t)/(x-t)v)dt. We say that f(x)is fuzzy Caputo fractional differentiable of order 0<v≤1atx0, if there exists an element(cDa+vf)(x0)∈Cℝℱ[a,b][a,b]such that for all0≤r≤1,h>0,
(26)(i)(cDa+vf)(x0)=limh→0+Φ(x0+h)⊖Φ(x0)h=limh→0+Φ(x0)⊖Φ(x0-h)h,
or
(27)(ii)(cDa+vf)(x0)=limh→0+Φ(x0)⊖Φ(x0+h)-h=limh→0+Φ(x0-h)⊖Φ(x0)-h,
or
(28)(iii)(cDa+vf)(x0)=limh→0+Φ(x0+h)⊖Φ(x0)h=limh→0+Φ(x0-h)⊖Φ(x0)-h,
or
(29)(iv)(cDa+vf)(x0)=limh→0+Φ(x0)⊖Φ(x0+h)-h=limh→0+Φ(x0)⊖Φ(x0-h)h.For the sake of simplicity, we say that the fuzzy-valued functionfisc[(1)-v]-differentiable if it is differentiable as in Definition 13 case (i), andfisc[(2)-v]-differentiable if it is differentiable as in Definition 13 case (ii), and so on for the other cases.Theorem 14 (see [63]).
Let0<v≤1andf∈ACℝℱ[a,b]; then the fuzzy Caputo fractional derivative exists almost everywhere on(a,b), and for all0≤r≤1one has
(30)(cDa+vf)(x;r)=[1Γ(1-v)∫axf-′r(t)dt(x-t)v,1Γ(1-v)∫ax‍f+′r(t)dt(x-t)v]=[(Ia+1-vDf-r)(x),(Ia+1-vDf+r)(x)],
whenfis (1)-differentiable, and
(31)(cDa+vf)(x;r)=[1Γ(1-v)∫ax‍f+′r(t)dt(x-t)v,1Γ(1-v)∫ax‍f-′r(t)dt(x-t)v,]=[(Ia+1-vDf+r)(x),(Ia+1-vDf-r)(x)],
when f is (2)-differentiable.Theorem 15 (fuzzy generalized Taylor’s formula [72]).
Letf(x)∈ACℝℱ[a,b](0,b], and suppose that cDkαf(x)∈Cℝℱ[a,b](0,b]fork=0,1,…,n+1where0<α<1,0≤x0≤xandx∈(0,b]. Then, one has
(32)[f(x)]r=[f_r(x),f¯r(x)],f_r(x)=∑i=0nxiαΓ(iα+1)cDiαf_r(0+)+cD(n+1)αfr(x0)Γ(nα+α+1)x(n+1)α,f¯r(x)=∑i=0nxiαΓ(iα+1)cDiαf¯r(0+)+cD(n+1)αfr(x0)Γ(nα+α+1)x(n+1)α,
where cDαf_r(0)=cDαf_r(x)|x=0,cDαf¯r(0)=cDαf¯r(x)|x=0.


## 2.2. Jacobi Polynomials

The well-known Jacobi polynomials associated with the parameters (α>-1,β>-1) (see, e.g., Luke [79] and Szegö [80]), are a sequence of polynomialsPi(α,β)(t)(i=0,1,…), each, respectively, of degreei. For using these polynomials on(0,L), we present the shifted Jacobi polynomials by implementing the change of variablet=(2x/L-1). Let the shifted Jacobi polynomials Pi(α,β)(2x/L-1)be denoted byPL,i(α,β)(x), satisfying the orthogonality relation
(33)∫0L‍PL,j(α,β)(x)PL,k(α,β)(x)wL(α,β)(x)dx=hk,
wherewL(α,β)(x)=xβ(L-x)αand
(34)hk={Lα+β+1Γ(k+α+1)Γ(k+β+1)(2k+α+β+1)k!Γ(k+α+β+1)i=j,0,i≠j.
The shifted Jacobi polynomialPL,i(α,β)(x)of degreeihas the form
(35)PL,i(α,β)(x)=∑k=0i‍(-1)i-kΓ(i+β+1)Γ(i+k+α+β+1)Γ(k+β+1)Γ(i+α+β+1)(i-k)!k!Lkxk,
where
(36)PL,i(α,β)(0)=(-1)iΓ(i+β+1)Γ(β+1)i!,PL,i(α,β)(L)=Γ(i+α+1)Γ(α+1)i!.
Also we can state the shifted Jacobi polynomial by the following recurrence relation:
(37)PL,i(α,β)(x)=(α+β+2i-1){(α2-β2+(2xL-1))×(α+β+2i)(α+β+2i-2)(2xL-1)}×(2i(α+β+i)(α+β+2i-2))-1×PL,i-1(α,β)(x)-(α+i-1)(β+i-1)(α+β+2i)i(α+β+i)(α+β+2i-2)×PL,i-2(α,β)(x)i=2,3,…,
wherePL,0(α,β)(x)=1andPL,1(α,β)(x)=((α+β+2)/2)(2x/L-1)+(α-β)/2.We notice that a functionu(x), square integrable in(0,L), can be expanded in terms of shifted Jacobi polynomials as
(38)u(x)=∑j=0+∞‍ajPL,j(α,β)(x),
where the coefficientsajare
(39)aj=1hj∫0L‍PL,j(α,β)(x)u(x)wL(α,β)(x)dx,j=0,1,…
For a given particular problem, only the first(N+1)-terms shifted Jacobi polynomials are used. Therefore, we have
(40)uN(x)≃∑j=0N‍ajPL,j(α,β)(x).Theorem 16 (see [59]).
The Caputo fractional derivative of ordervof the shifted Jacobi polynomials of degreeiis obtained from
(41)DvPL,i(α,β)(x)=∑j=0∞‍Sv(i,j,α,β)PL,i(α,β)(x),i=⌈v⌉,⌈v⌉+1,…,
where
(42)Sv(i,j,α,β)=∑k=⌈v⌉i(-1)i-kLα+β-v+1Γ(j+β+1)×Γ(i+β+1)Γ(i+k+α+β+1)×(hjΓ(j+k+α+β+1)Γ(k+β+1)×Γ(i+α+β+1)Γ(k-v+1)(i-k)!)-1×∑l=0j(-1)j-lΓ(j+l+α+β+1)Γ(α+1)×Γ(l+k+β-v+1)×(Γ(l+β+1)Γ(l+k+α+β-v+2)(j-l)!l!)-1.


## 3. Pharmacokinetics Model Equation

### 3.1. Drug Assimilation into the Blood

The drug dissolves in the gastrointestinal tract (GI), and each ingredient is diffused into the bloodstream. They are carried to the locations in which they act and are removed from the blood by the kidneys and liver. Generally, the problem of drug assimilation into the body can be considered as a two-compartment model: GI-tract and the bloodstream [16]. Different compartments and the input/output of the model are depicted in Figure 1.Figure 1
Schematic of input/output compartment for drug assimilation.For each compartment by applying the balance law we can obtain(43){rateofchangeofdruginGItract}={rateofdrugintake}-{ratedrugleavesGI-tract}{rateofchangeofdruginblood}={ratedrugentersblood}-{ratedrugleavesblood}.
In this study, we consider a case of a single cold pill. Also, there is no ingestion of the drug except that which occurs initially.


### 3.2. Case of a Single Cold Pill

Let us considerx(t)to be the amount of drug in the GI-tract at timetandy(t), the amount in the bloodstream at timet. In the GI-tract we suppose that the pill is to be swallowed, and so after this event (over subsequent time) we have nothing more entering the GI-tract. The pill dissolves and diffuses into the bloodstream from the GI-tract. So, the GI-tract is only an output term. Assuming that the output rate is proportional to the amount of drug in the bloodstream, which is in fact proportional to the GI-tract drug concentration, then [16],
(44)dydt=-k1x,x(0)=x0,
wherex0is the amount of a drug in the pill andk1is a positive coefficient. We suppose thaty(0)=0, which means that the initial amount of the drug in the bloodstream is zero. As the drug diffuses from the GI-tract, the level increases, and as the kidneys and liver remove it, the level of drug gradually decreases. Thus,
(45)dydt=k1x-k2y,y(0)=0,
withk2another positive constant. Decongestant and an antihistamine are the component of the cold pill, and the coefficient of proportionality (k1andk2) is different for the different component drugs in the pill.In order to obtain the growth and decay of antihistamine levels in the GI tract and bloodstream, rearrange the first rate equation, multiply by the integrating factorek1t, integrate, and then use the initial data to obtain
(46)X=Ae(-k1t).
Insert this formula forx(t)into the second rate equation, which becomes
(47)dy(t)dt+k2y(t)=k1Ae-k1ty(0)=0.
In this paper, we are concerned with fractional time derivatives regarding (17) and (18). Having defined cDv, we can now turn to fractional differential equations and systems of fractional differential equations, which will be used to specify PKPD models and will need to be solved over an interval[0,t], in accordance with appropriate initial conditions.A typical feature of differential equations (both classical and fractional) is the need to specify additional conditions in order to produce a unique solution. For the case of Caputo fractional differential equations, these additional conditions are just the initial conditions which are similarly required by classical ODEs. In contrast, for Riemann-Liouville fractional differential equations, these additional conditions constitute certain fractional derivatives (and/or integrals) of the unknown solution at the initial pointt=0[81], which are functions oft. These initial conditions are not physical; furthermore, it is not clear how such quantities are to be measured from experiment, say, so that they can be appropriately assigned in an analysis [82]. If for no other reason, the need to solve fractional differential equations is justification enough for choosing Caputo’s definition for fractional differentiation over the more commonly used (at least in mathematical analysis) definition of Liouville and Riemann, and this is the operator that we choose to use in the following.We start by representing drug concentration in the effect compartment by the (Caputo) fractional differential equation:(48)cDvy(t)+k2y(t)=k1Ae-k1ty(0)=0.In the standard direct action model, the effect at timet,Y(t), is expressed by an arbitrary (memory-less) function of drug concentration in the effect site at timet,G(y(t)); however to generate a wider class of relationships, we assume that the effect at timet∈[0,1]is related to the fuzzy Caputo fractional derivative ofy(t). So, we have
(49)cDvy(t)+k2y(t)=k1Ae-k1ty(0;r)=[y_0r,y¯0r],
in whichy(x):Lℝℱ[0,1]∩Cℝℱ[0,1]is a continuous fuzzy-valued function andcD0+vdenotes the fuzzy Caputo fractional derivative of orderv∈[0,1].Remark 17.
In this paper, the drug concentration is modeled by an oscillation-relaxation fuzzy fractional differential equation. So the right hand side of (49) is determined based on the model parameters. It should mention here that the proposed method can be extended easily for solving other types of FFDEs with more complicated right hand side functions. For more details, one can refer to (33) in [65] and [57–59].We have assumed thatk1≠k2, an assumption that is justified by the pharmaceutical data. For the “average” person, a pharmaceutical company estimates that the values of the rate constants for the antihistamine in the cold pills arek1=0.6931(hour)-1andk2=0.0231(hour)-1. It can be observed from (46) that level of antihistamine in the bloodstream increases as the time increases and saturate in a maximum value of antihistamine; however, (47) can conclude that as the time increases the amount of antihistamine in GI-tract declines and reaches a minimum value. We now consider two different cases of patients who are not average. We often definek2as clearance coefficient of medication from the blood. The value ofk2often varies from old and sick patients than young and healthy cases. This means that the level of medication in the blood may become and then remain excessively high with a standard dosage for the cases thatk2is much lower, which is normally observed in old and sick cases.We investigate sensitivity of the medication over a 24-hour period by keeping the value ofk1fixed at some value such as 1.386 but settingk2=0.01386,0.06386,0.1386,0.6386,1.386(A=1). Different values ofk2correspond to people of different ages and states of health. For all cases as it is expected that, for the first few hours, the amount of antihistamine in the bloodstream increases and then declines gradually. However, for the cases with lower value of clearance coefficient (old and sick), the maximum value of the medication in bloodstream was much higher than the cases with highk2, and it did not decline for remaining times. This means that the level of medication in the bloodstream stays at high level for a long time, and it could not be absorbed from the blood. The same trend is observed for the case thatk2is constant (k2=0.0231) and the value ofk1varies (k1=0.06931,0.11,0.3,0.6931,1.0, and1.5). In this case, thek1is a constant for the GI tract which is analogous to thek2for the bloodstream.


## 3.1. Drug Assimilation into the Blood

The drug dissolves in the gastrointestinal tract (GI), and each ingredient is diffused into the bloodstream. They are carried to the locations in which they act and are removed from the blood by the kidneys and liver. Generally, the problem of drug assimilation into the body can be considered as a two-compartment model: GI-tract and the bloodstream [16]. Different compartments and the input/output of the model are depicted in Figure 1.Figure 1
Schematic of input/output compartment for drug assimilation.For each compartment by applying the balance law we can obtain(43){rateofchangeofdruginGItract}={rateofdrugintake}-{ratedrugleavesGI-tract}{rateofchangeofdruginblood}={ratedrugentersblood}-{ratedrugleavesblood}.
In this study, we consider a case of a single cold pill. Also, there is no ingestion of the drug except that which occurs initially.


## 3.2. Case of a Single Cold Pill

Let us considerx(t)to be the amount of drug in the GI-tract at timetandy(t), the amount in the bloodstream at timet. In the GI-tract we suppose that the pill is to be swallowed, and so after this event (over subsequent time) we have nothing more entering the GI-tract. The pill dissolves and diffuses into the bloodstream from the GI-tract. So, the GI-tract is only an output term. Assuming that the output rate is proportional to the amount of drug in the bloodstream, which is in fact proportional to the GI-tract drug concentration, then [16],
(44)dydt=-k1x,x(0)=x0,
wherex0is the amount of a drug in the pill andk1is a positive coefficient. We suppose thaty(0)=0, which means that the initial amount of the drug in the bloodstream is zero. As the drug diffuses from the GI-tract, the level increases, and as the kidneys and liver remove it, the level of drug gradually decreases. Thus,
(45)dydt=k1x-k2y,y(0)=0,
withk2another positive constant. Decongestant and an antihistamine are the component of the cold pill, and the coefficient of proportionality (k1andk2) is different for the different component drugs in the pill.In order to obtain the growth and decay of antihistamine levels in the GI tract and bloodstream, rearrange the first rate equation, multiply by the integrating factorek1t, integrate, and then use the initial data to obtain
(46)X=Ae(-k1t).
Insert this formula forx(t)into the second rate equation, which becomes
(47)dy(t)dt+k2y(t)=k1Ae-k1ty(0)=0.
In this paper, we are concerned with fractional time derivatives regarding (17) and (18). Having defined cDv, we can now turn to fractional differential equations and systems of fractional differential equations, which will be used to specify PKPD models and will need to be solved over an interval[0,t], in accordance with appropriate initial conditions.A typical feature of differential equations (both classical and fractional) is the need to specify additional conditions in order to produce a unique solution. For the case of Caputo fractional differential equations, these additional conditions are just the initial conditions which are similarly required by classical ODEs. In contrast, for Riemann-Liouville fractional differential equations, these additional conditions constitute certain fractional derivatives (and/or integrals) of the unknown solution at the initial pointt=0[81], which are functions oft. These initial conditions are not physical; furthermore, it is not clear how such quantities are to be measured from experiment, say, so that they can be appropriately assigned in an analysis [82]. If for no other reason, the need to solve fractional differential equations is justification enough for choosing Caputo’s definition for fractional differentiation over the more commonly used (at least in mathematical analysis) definition of Liouville and Riemann, and this is the operator that we choose to use in the following.We start by representing drug concentration in the effect compartment by the (Caputo) fractional differential equation:(48)cDvy(t)+k2y(t)=k1Ae-k1ty(0)=0.In the standard direct action model, the effect at timet,Y(t), is expressed by an arbitrary (memory-less) function of drug concentration in the effect site at timet,G(y(t)); however to generate a wider class of relationships, we assume that the effect at timet∈[0,1]is related to the fuzzy Caputo fractional derivative ofy(t). So, we have
(49)cDvy(t)+k2y(t)=k1Ae-k1ty(0;r)=[y_0r,y¯0r],
in whichy(x):Lℝℱ[0,1]∩Cℝℱ[0,1]is a continuous fuzzy-valued function andcD0+vdenotes the fuzzy Caputo fractional derivative of orderv∈[0,1].Remark 17.
In this paper, the drug concentration is modeled by an oscillation-relaxation fuzzy fractional differential equation. So the right hand side of (49) is determined based on the model parameters. It should mention here that the proposed method can be extended easily for solving other types of FFDEs with more complicated right hand side functions. For more details, one can refer to (33) in [65] and [57–59].We have assumed thatk1≠k2, an assumption that is justified by the pharmaceutical data. For the “average” person, a pharmaceutical company estimates that the values of the rate constants for the antihistamine in the cold pills arek1=0.6931(hour)-1andk2=0.0231(hour)-1. It can be observed from (46) that level of antihistamine in the bloodstream increases as the time increases and saturate in a maximum value of antihistamine; however, (47) can conclude that as the time increases the amount of antihistamine in GI-tract declines and reaches a minimum value. We now consider two different cases of patients who are not average. We often definek2as clearance coefficient of medication from the blood. The value ofk2often varies from old and sick patients than young and healthy cases. This means that the level of medication in the blood may become and then remain excessively high with a standard dosage for the cases thatk2is much lower, which is normally observed in old and sick cases.We investigate sensitivity of the medication over a 24-hour period by keeping the value ofk1fixed at some value such as 1.386 but settingk2=0.01386,0.06386,0.1386,0.6386,1.386(A=1). Different values ofk2correspond to people of different ages and states of health. For all cases as it is expected that, for the first few hours, the amount of antihistamine in the bloodstream increases and then declines gradually. However, for the cases with lower value of clearance coefficient (old and sick), the maximum value of the medication in bloodstream was much higher than the cases with highk2, and it did not decline for remaining times. This means that the level of medication in the bloodstream stays at high level for a long time, and it could not be absorbed from the blood. The same trend is observed for the case thatk2is constant (k2=0.0231) and the value ofk1varies (k1=0.06931,0.11,0.3,0.6931,1.0, and1.5). In this case, thek1is a constant for the GI tract which is analogous to thek2for the bloodstream.


## 4. Description of the Method

In this study, by developing the Jacobi polynomial approximation [57, 58, 65] with the help of the matrix operations, the tau method, and the fuzzy Caputo fractional derivative, we obtain an approximate solution of the problem (49) for different values ofk1andk2. As it is known, the existence and uniqueness of the FFDEs have been introduced in [60, 61, 63]. Here, we assume that the shifted Jacobi polynomials are defined on(0,1), so we have [65]
(50)u(x)=∑i=0+∞‍*ai⊙P1,i(α,β)(x),
where the fuzzy coefficientsaiare gained by
(51)ai=1hi∫01‍P1,i(α,β)(x)⊙u(x)⊙w1(α,β)(x)dx,i=0,1,…,
andw1(α,β)(x)=(1-x)α⊙xβ,u∈Lpℝℱ[0,1]∩Cℝℱ[0,1], andP1,i(α,β)(x)is as the same as the shifted Jacobi polynomials presented in Section 2.2, and∑‍*means addition with respect to⊕inℝℱ.Remark 18.
In the remainder of paper, for more simplicity, we considerPi(α,β)(x)instead ofP1,i(α,β)(x)andw(α,β)(x)=(1-x)α⊙xβinstead ofw1(α,β)(x)=(1-x)α⊙xβ.Remark 19.
Practically, only the first(N+1)-terms shifted Jacobi polynomials are taken into consideration. So, we have
(52)u(x)≃uN(x)=∑i=0N‍a*i⊙Pi(α,β)(x)=aT⊙ΦN(x),
where the fuzzy shifted Jacobi coefficient vectorATand shifted Jacobi polynomials vectorΦN+1(x)are presented by
(53)aT=[a0,a1,…,aN],ΦN(x)=[P0(α,β)(x),P1(α,β)(x),…,PN(α,β)(x)]T.
We can state the fuzzy approximate function (52) in the parametric form as follows.Definition 20 (see [65]).
Letu(x)∈Lpℝℱ[0,1]∩Cℝℱ[0,1]; the approximation of fuzzy-valued functionu(x)in the parametric form is
(54)ur(x)≃uNr(x)=[∑i=0N‍ai,-rPi(α,β)(x),∑i=0N‍ai,+rPi(α,β)(x)],0≤r≤1.We call thatSpan{Pn(α,β)(x):n≥0}forms a complete orthogonal system inL2,w(α,β)ℝℱ(0,1). Hence, we set
(55)𝕊N(0,1)=Span{P0(α,β)(x),P1(α,β)(x),…,PN(α,β)(x)}.
The shifted Jacobi tau method to (49) is to obtainuN∈SN(0,1)such that
(56)(DvuN,Pk(α,β)(x))w(α,β)+k2(uN,Pk(α,β)(x))w(α,β)=(k1Ae-k1t,Pk(α,β)(x))w(α,β),k=0,1,…,N-1,uN(0)=d0,
wherew(α,β)(x)=(1-x)α⊙xβand(u,v)w(α,β)=∫01‍w(α,β)(x)⊙u(x)⊙v(x)dxdenotes the fuzzy inner product inL2,w(α,β)ℝℱ(0,1)andd0is the fuzzy initial value of the problem.Let us define the following notations:(57)fk=(k1Ae-k1t,Pk(α,β)(x))w(α,β)k=0,1,…,N-1,f=(f0,f1,…,fN-1,d0)T.
Hence, the variational formulation of (56) according to Relation (14) in [65], by means of a typical tau method like in the crisp context [29] and (52), is equivalent to
(58)∑j=0N‍aj⊙[(D(v)Pj(α,β)(x),Pk(α,β)(x))w(α,β)+k2(Pj(α,β)(x),Pk(α,β)(x))w(α,β)]=(k1Ae-k1t,Pk(α,β)(x))w(α,β),k=0,1,…,N-1,∑j=0N‍aj⊙Pj(α,β)(0)=d0.
Denoting
(59)A=(akj)0<k,j<N,C=(ckj)0<k,j<N,
we investigate that (58) is equivalent to the matrix system
(60)(A+k2C)⊙a=f.
The elements of the matrices mentioned above are determined in the next theorem.Theorem 21 (see [59]).
Let us denoteakj=(D(v)Pj(α,β)(x),Pk(α,β)(x))w(α,β)(0⩽k⩽N-1,0⩽j⩽N),akj=Dk-NPj(α,β)(0)(k=N,0⩽j⩽N) and ckj=(Pj(α,β)(x),Pk(α,β)(x))w(α,β)(0⩽k⩽N-1,0⩽j⩽N); then the nonzero elements ofakjandckjare given as
(61)akj=hkSv(j,k,α,β),0⩽k⩽N-1,1⩽j⩽N,akj=(-1)j-k+NΓ(j+β+1)(j+α+β+1)k-NΓ(j-k+N+1)Γ(k-N+1+β),k=N,0⩽j⩽N,ckj=hk,0⩽k=j⩽N-1.


### 4.1. Error Analysis

In this section, error analysis of the method will be presented for the FFDEs. Firstly, an upper bound of the absolute errors will be given for the technique by using generalized Taylor formula. Secondly, an error bound will be introduced for the approximation of the fuzzy fractional Caputo’s derivative using Jacobi polynomials.Lemma 22.
LetyN,v(t)andf(t)be the Jacobi approximate solution (40) and the exact solution of (48), respectively. IfD0+kvf(t)∈C(0,b],k=0,1,…,N+1, then
(62)|f(t)-yN,v(t)|⩽|RNv(t,0)|+|fTNi(t)-yN,v(t)|,wherefTNiandRNvare defined according to Theorem 10.Proof.
SinceD0+kvf(t)∈C(0,b],k=0,1,…,N+1,fcan be expanded to the fractional Taylor series
(63)fTNi(t)=∑i=0N‍tivΓ(iv+1)[D0+ivf(t)]t=0+,
and its reminder term is
(64)RNv(t,0)=t(N+1)vΓ((N+1)v+1)[D0+(N+1)vf(t)]t=ξ,0⩽ξ⩽t,∀t∈(0,b].
Therefore,
(65)f(t)-fTNi(t)=RNv(t,0).
From (65) and triangle inequality, we obtain
(66)|f(t)-yN,v(t)|=|f(t)-yN,v(t)+fTNi(t)-fTNi(t)|⩽|f(t)-fTNi(t)|+|fTNi(t)-yN,v(t)|=|RNv(t,0)|+|fTNi(t)-yN,v(t)|.
Therefore, an upper bound of the absolute errors is obtained for the method in the crisp cases.Now, we provide an upper bound for the absolute errors of the fuzzy approximate solution by using the proposed method.Theorem 23.
Letu(t)∈Lpℝℱ[0,1]∩Cℝℱ[0,1]be the fuzzy exact solution of (49), anduN(t)is the best fuzzy Jacobi approximate function (52), and suppose thatcDkvu(t)∈Cℝℱ[0,b]fork=0,1,…,N+1where0<v<1andt∈[0,b]. Then, one has
(67)[u(t)]r=[u_r(t),u¯r(t)],[uN(t)]r=[u_Nr(t),u¯Nr(t)],|u_(t)-u_N(t)|⩽|R_Nv(t,0)|+|u_TNi(t)-u_N(t)|,|u¯(t)-u¯N(t)|⩽|R¯Nv(t,0)|+|u¯TNi(t)-u¯N(t)|,
where regarding Theorem 15, one has
(68)[fTNi(t)]r=[f_TNi(t;r),f¯TNi(t;r)]=[∑i=0N‍tivΓ(iv+1)[D0+ivf_(t)]t=0+,∑i=0N‍tivΓ(iv+1)[D0+ivf¯(t)]t=0+],[RNv(t,0)]r=[R_Nv,r(t,0),R¯Nv,r(t,0)]=[cD(N+1)vf_r(0+)Γ(Nv+v+1)t(N+1)v,cD(N+1)vf¯r(0+)Γ(Nv+v+1)t(N+1)vcD(N+1)vf_r(0+)Γ(Nv+v+1)].Proof.
It is straightforward from Lemma22 and the fuzzy function definition.Also in the following theorem, according to the Relation (14) in [65], an upper bound for the absolute errors of the approximate function of fuzzy fractional Caputo’s derivative is provided.Theorem 24 (see [65]).
Assume that the error function of fuzzy Caputo fractional derivative operator by using the shifted Jacobi polynomials,Ek,v, is continuously fuzzy differentiable for0<x0≤x,x∈(0,1]. Additionally,Ek,v∈Cℝℱ[x0,1]and0<v<1; then the error bound is given by
(69)D*(ℝℱDvΦ(x),ℝℱD(v)Φ(x))≤x0-v(2)!|Γ(1-v)|(S)2(β+22)B(α+1,β+1).


## 4.1. Error Analysis

In this section, error analysis of the method will be presented for the FFDEs. Firstly, an upper bound of the absolute errors will be given for the technique by using generalized Taylor formula. Secondly, an error bound will be introduced for the approximation of the fuzzy fractional Caputo’s derivative using Jacobi polynomials.Lemma 22.
LetyN,v(t)andf(t)be the Jacobi approximate solution (40) and the exact solution of (48), respectively. IfD0+kvf(t)∈C(0,b],k=0,1,…,N+1, then
(62)|f(t)-yN,v(t)|⩽|RNv(t,0)|+|fTNi(t)-yN,v(t)|,wherefTNiandRNvare defined according to Theorem 10.Proof.
SinceD0+kvf(t)∈C(0,b],k=0,1,…,N+1,fcan be expanded to the fractional Taylor series
(63)fTNi(t)=∑i=0N‍tivΓ(iv+1)[D0+ivf(t)]t=0+,
and its reminder term is
(64)RNv(t,0)=t(N+1)vΓ((N+1)v+1)[D0+(N+1)vf(t)]t=ξ,0⩽ξ⩽t,∀t∈(0,b].
Therefore,
(65)f(t)-fTNi(t)=RNv(t,0).
From (65) and triangle inequality, we obtain
(66)|f(t)-yN,v(t)|=|f(t)-yN,v(t)+fTNi(t)-fTNi(t)|⩽|f(t)-fTNi(t)|+|fTNi(t)-yN,v(t)|=|RNv(t,0)|+|fTNi(t)-yN,v(t)|.
Therefore, an upper bound of the absolute errors is obtained for the method in the crisp cases.Now, we provide an upper bound for the absolute errors of the fuzzy approximate solution by using the proposed method.Theorem 23.
Letu(t)∈Lpℝℱ[0,1]∩Cℝℱ[0,1]be the fuzzy exact solution of (49), anduN(t)is the best fuzzy Jacobi approximate function (52), and suppose thatcDkvu(t)∈Cℝℱ[0,b]fork=0,1,…,N+1where0<v<1andt∈[0,b]. Then, one has
(67)[u(t)]r=[u_r(t),u¯r(t)],[uN(t)]r=[u_Nr(t),u¯Nr(t)],|u_(t)-u_N(t)|⩽|R_Nv(t,0)|+|u_TNi(t)-u_N(t)|,|u¯(t)-u¯N(t)|⩽|R¯Nv(t,0)|+|u¯TNi(t)-u¯N(t)|,
where regarding Theorem 15, one has
(68)[fTNi(t)]r=[f_TNi(t;r),f¯TNi(t;r)]=[∑i=0N‍tivΓ(iv+1)[D0+ivf_(t)]t=0+,∑i=0N‍tivΓ(iv+1)[D0+ivf¯(t)]t=0+],[RNv(t,0)]r=[R_Nv,r(t,0),R¯Nv,r(t,0)]=[cD(N+1)vf_r(0+)Γ(Nv+v+1)t(N+1)v,cD(N+1)vf¯r(0+)Γ(Nv+v+1)t(N+1)vcD(N+1)vf_r(0+)Γ(Nv+v+1)].Proof.
It is straightforward from Lemma22 and the fuzzy function definition.Also in the following theorem, according to the Relation (14) in [65], an upper bound for the absolute errors of the approximate function of fuzzy fractional Caputo’s derivative is provided.Theorem 24 (see [65]).
Assume that the error function of fuzzy Caputo fractional derivative operator by using the shifted Jacobi polynomials,Ek,v, is continuously fuzzy differentiable for0<x0≤x,x∈(0,1]. Additionally,Ek,v∈Cℝℱ[x0,1]and0<v<1; then the error bound is given by
(69)D*(ℝℱDvΦ(x),ℝℱD(v)Φ(x))≤x0-v(2)!|Γ(1-v)|(S)2(β+22)B(α+1,β+1).


## 5. Numerical Results

In this section, the fuzzy fractional PKPD model (49) is solved for different values ofk1andk2by using the JT method presented in Section 4. We have performed all numerical computations with a computer program written in MATLAB. Also, absolute errors between fuzzy approximate solution[yN]r=[y_Nr,y¯Nr]and the corresponding exact solutions y(t;r)=[y_(t;r),y¯(t;r)], that is, [Ne]r=[|y_Nr-y_r|,|y¯Nr-y¯r|], are considered.Now, we recall the FFDE (49) as follows:
(70)cDvy(t)+k2y(t)=k1Ae-k1ty(0;r)=[-1+r,1-r],
in whichy(x):Lℝℱ[0,1]∩Cℝℱ[0,1]is a continuous fuzzy-valued function andcD0+vdenotes the fuzzy Caputo fractional derivative of orderv∈[0,1].Two cases are considered for solving by the proposed technique. In the first case, we assume thatk2is unchangeable and different values ofk1are substituted in (49) to get the fuzzy approximate solution. Conversely, in the second case we try to solve the problem such thatk1has an invariable value andk2is varied. The problem is analyzed for both cases in details.Case I.
Let us consider from Section3 thatk2=0.0231, but letk1vary (e.g.,0.6931,0.11, and0.3). So (49), with the assumption, namely,k2=0.0231and DoseA=1, is as follows:
(71)cDvy(t)+0.0231y(t)=0.6931e-0.6931t,y(0;r)=[-1+r,1-r].By using Theorem7 for the above equation under c[(1)-v] differentiability, we have the following systems:
(72)cDvy_(t;r)+0.0231y_(t;r)=0.6931e-0.6931t,0<v≤1,y_(0;r)=-1+r,0<r≤1,cDvy¯(t;r)+0.0231y¯(t;r)=0.6931e-0.6931t,0<v≤1,y¯(0;r)=1-r,0<r≤1.Solving (72) leads to determining the exact solution of (71) as follows:
(73)Y_(t;r)=(-1+r)Ev,1[-0.0231tv]+∫0t‍(t-x)v-1Ev,v[-0.0231(t-x)v]×(0.6931e-0.6931x)dx,0<v≤1,Y¯(t;r)=(1-r)Ev,1[-0.0231tv]+∫0t‍(t-x)v-1Ev,v[-0.0231(t-x)v]×(0.6931e-0.6931x)dx,0<r≤1.We seek the fuzzy approximate solutions by applying the technique described in Section4 withN=2and as
(74)y2(t)=a0⊙P0(α,β)(t)+a1⊙P1(α,β)(t)+a2⊙P0(α,β)(t).
Here, we have
(75)D0.85=(0001.86390.3901-0.1755-0.39014.52670.8696),
whereα=β=0. Regarding (54), we have
(76)y_2r(t)=a_0r+a_1r(2x-1)+a_2r(6x2-6x+1)y¯2r(t)=a¯0r+a¯1r(2x-1)+a¯2r(6x2-6x+1).By solving the fuzzy linear system (60) in the parametric form, one can get fuzzy unknown coefficients in (76) as
(77)a_0=-0.5757,a_1=0.2624,a_2=-0.0619,a¯0=1.2008,a¯1=0.2408,a¯2=-0.0600,
where we assume thatr-cut=0.1in (77).The comparison between absolute errors of differentk1obtained by our method is shown in Table 1. In Figure 2 logarithmic plot of absolute error with different values ofαandβwas obtained, and different number of Jacobi functions was experienced for solving Case I by using the proposed method in Figure 3. Finally, the approximate fuzzy solutions are shown for different values ofk1with the fractional order v=0.85 in Figure 4.Table 1
The results of the proposed method for Case I withv=0.85, α=β=0, and N=8.
r
k
1
=
0.6931
k
1
=
0.11
k
1
=
0.3
k
1
=
0.6931
k
1
=
0.11
k
1
=
0.3
N
_
e
r
N
_
e
r
N
_
e
r
N
¯
e
r
N
¯
e
r
N
¯
e
r
0
9.0140
e
-
5
1.6623
e
-
5
4.0452
e
-
5
8.4377
e
-
5
1.0860
e
-
5
3.4689
e
-
5
0.1
8.9852
e
-
5
1.6335
e
-
5
4.0164
e
-
5
8.4665
e
-
5
1.1148
e
-
5
3.4977
e
-
5
0.2
8.9564
e
-
5
1.6046
e
-
5
3.9876
e
-
5
8.4953
e
-
5
1.1436
e
-
5
3.5265
e
-
5
0.3
8.9276
e
-
5
1.5758
e
-
5
3.9588
e
-
5
8.5242
e
-
5
1.1724
e
-
5
3.5553
e
-
5
0.4
8.8988
e
-
5
1.5470
e
-
5
3.9300
e
-
5
8.5530
e
-
5
1.2012
e
-
5
3.5842
e
-
5
0.5
8.8700
e
-
5
1.5182
e
-
5
3.9011
e
-
5
8.5818
e
-
5
1.2300
e
-
5
3.6130
e
-
5
0.6
8.8411
e
-
5
1.4894
e
-
5
3.8723
e
-
5
8.6106
e
-
5
1.2588
e
-
5
3.6418
e
-
5
0.7
8.8123
e
-
5
1.4606
e
-
5
3.8435
e
-
5
8.6394
e
-
5
1.2877
e
-
5
3.6706
e
-
5
0.8
8.7835
e
-
5
1.4317
e
-
5
3.8147
e
-
5
8.6682
e
-
5
1.3165
e
-
5
3.6994
e
-
5
0.9
8.7547
e
-
5
1.4029
e
-
5
3.7859
e
-
5
8.6971
e
-
5
1.3453
e
-
5
3.7282
e
-
5
1
8.7259
e
-
5
1.3741
e
-
5
3.7571
e
-
5
8.7259
e
-
5
1.3741
e
-
5
3.7571
e
-
5Figure 2
The absolute errors for differentαandβwithN=9,k1=0.3andv=0.75, Case I.Figure 3
The absolute errors for different valuesNwithv=0.95,k1=0.11andα=β=0, Case I.Figure 4
The fuzzy approximate solution of Case I, for different value ofk1,α=β=0.5,N=8.The absolute error value forr-cut varied from 0 to 1 for different value ofk1which is calculated in Table 1. As it can be observed at a constantr-cut by increasing the value ofk1which is the ratio of variation of the drug in the bloodstream to the amount of drug in the GI-tract the value of absolute error increases. This is analogous to the previous reports for various values of thek1[12]. It should be noted here that the value ofNeris directly proportional to the variation of exact solutions and corresponding fuzzy approximate solution. As a result, increase ofNerwith increasingk1is also expected.Case II.
In this case, we assume thatk1=1.386and is not variable, but k2 takes different values (i.e.,k2=0.01386,0.06386, and 0.1386) and DoseA=1. So (70) alters with, forexample,k2=0.01386as follows:
(78)cDvy(t)+0.01386y(t)=1.386e-1.386t,y(0;r)=[-1+r,1-r],
with the exact solution as:
(79)Y_(t;r)=(-1+r)Ev,1[-0.01386tv]+∫0t‍(t-x)v-1Ev,v[-0.01386(t-x)v]×(1.386e-1.386x)dx,0<v≤1,Y¯(t;r)=(1-r)Ev,1[-0.01386tv]+∫0t‍(t-x)v-1Ev,v[-0.01386(t-x)v]×(1.386e-1.386x)dx,0<r≤1.By applying the JT method in Section4, we can get the fuzzy unknown coefficients{aj}j=0Nby using (58) which is in the parametric form as follows:
(80)∑j=0Na_j[(D(v)Pj(α,β)(x),Pk(α,β)(x))w(α,β)+0.01386(Pj(α,β)(x),Pk(α,β)(x))w(α,β)]=(1.386Ae-1.386t,Pk(α,β)(x))w(α,β),k=0,1,…,N-1,∑j=0N‍a_jPj(α,β)(0)=-1+r,∑j=0Na¯j[(D(v)Pj(α,β)(x),Pk(α,β)(x))w(α,β)+0.01386(Pj(α,β)(x),Pk(α,β)(x))w(α,β)]=(1.386Ae-1.386t,Pk(α,β)(x))w(α,β),k=0,1,…,N-1,∑j=0N‍a¯jPj(α,β)(0)=1-r,
whereD(v)is given as
(81)D0.75=(0002.69290.5524-0.1755-1.24294.22411.1048),withα=β=0.5andN=2. By substituting the above matric in (80) we can reach the fuzzy linear algebraic system (60) that can solve easily to determine the fuzzy coefficients as shown in Table 2.Table 2
The fuzzy coefficients for Case II withv=0.75, α=β=0.5, N=2, and k2=0.01386.
r
a
_
0
a
_
1
a
_
2
a
¯
0
a
¯
1
a
¯
2
0
−0.42793
0.24218
−0.11135
1.55482
0.23234
−0.11003
0.1
−0.32879
0.24169
−0.11128
1.45568
0.23283
−0.11009
0.2
−0.22965
0.24120
−0.11121
1.35654
0.23332
−0.11016
0.3
−0.13052
0.24071
−0.11115
1.25741
0.23382
−0.11022
0.4
−0.03138
0.24021
−0.11108
1.15827
0.23431
−0.11029
0.5
0.06775
0.23972
−0.11102
1.05913
0.23480
−0.11036
0.6
0.16689
0.23923
−0.11095
0.95999
0.23529
−0.11042
0.7
0.26603
0.23874
−0.11088
0.86085
0.23578
−0.11049
0.8
0.36516
0.23824
−0.11082
0.76172
0.23628
−0.11055
0.9
0.46430
0.23775
−0.11075
0.66258
0.23677
−0.11062
1
0.56344
0.23726
−0.11069
0.56344
0.23726
−0.11069We compared JT method results for different values ofk2and the outcomes are tabulated in Table 3. The absolute errors of the proposed method for this case are exhibited in Figure 5 with four choices ofαand β. Clearly, the best approximation is achieved when the values ofα=0andβ=0. Additionally, we compared the approximate solutions obtained by the present method atN=4,7,9, and11which is shown in Figure 6. From Figure 6, one can conclude that with increasing the number of Jacobi polynomials, the absolute errors are decreasing dramatically. Finally, Figure 7 shows the numerical results for different values ofk2atv=0.75.It can be seen that for all the values ofk2the approximate solution is the fuzzy number.Table 3
The results of the proposed method for Case II withv=0.95, α=β=0.5, and N=7.
r
k
2
=
0.01386
k
2
=
0.06386
k
2
=
0.1386
k
2
=
0.01386
k
2
=
0.6386
k
2
=
0.1386
N
_
e
r
N
_
e
r
N
_
e
r
N
¯
e
r
N
¯
e
r
N
¯
e
r
0
5.2740
e
-
4
5.2151
e
-
4
5.1213
e
-
4
5.1711
e
-
4
4.7626
e
-
4
4.2036
e
-
4
0.1
5.2688
e
-
4
5.1925
e
-
4
5.0754
e
-
4
5.1763
e
-
4
4.7852
e
-
4
4.2495
e
-
4
0.2
5.2637
e
-
4
5.1698
e
-
4
5.0295
e
-
4
5.1814
e
-
4
4.8078
e
-
4
4.2953
e
-
4
0.3
5.2585
e
-
4
5.1472
e
-
4
4.9836
e
-
4
5.1866
e
-
4
4.8304
e
-
4
4.3412
e
-
4
0.4
5.2534
e
-
4
5.1246
e
-
4
4.9377
e
-
4
5.1917
e
-
4
4.8531
e
-
4
4.3871
e
-
4
0.5
5.2482
e
-
4
5.1019
e
-
4
4.8919
e
-
4
5.1968
e
-
4
4.8757
e
-
4
4.4330
e
-
4
0.6
5.2431
e
-
4
5.0793
e
-
4
4.8460
e
-
4
5.2020
e
-
4
4.8983
e
-
4
4.4789
e
-
4
0.7
5.2380
e
-
4
5.0567
e
-
4
4.8001
e
-
4
5.2071
e
-
4
4.9209
e
-
4
4.5248
e
-
4
0.8
5.2328
e
-
4
5.0341
e
-
4
4.7542
e
-
4
5.2123
e
-
4
4.9436
e
-
4
4.5707
e
-
4
0.9
5.2277
e
-
4
5.0114
e
-
4
4.7083
e
-
4
5.2174
e
-
4
4.9662
e
-
4
4.6165
e
-
4
1
5.2225
e
-
4
4.9888
e
-
4
4.6624
e
-
4
5.2225
e
-
4
4.9888
e
-
4
4.6624
e
-
4Figure 5
The absolute errors for differentαandβwithN=8,k2=0.1386, andv=0.95, Case II.Figure 6
The absolute errors for different valuesNwithv=0.85,k2=0.1386, andα=β=0, Case II.Figure 7
The fuzzy approximate solution of Case II, for different value ofk2,α=0,β=0.5, andN=9.In Table3, the value of absolute error forr-cut varied from 0 to 1 for different value ofk2is calculated. At a constantr-cut value, a descending trend is observed with increasing the value ofk2. This trend is opposed to the trend observed in Table 1 for different value ofk1.k2is clearance constant and it is a coefficient fory(t)which is the amount of drug in the bloodstream in (44). We already explained that by decreasing the value ofk2,which happens in the case of old and sick patient, the drug absorption from the bloodstream by kidney and liver also decreases and the drug stays in the bloodstream for longer time. Since the constant value ofk2appears with negative sign in the main differential equation as we expected, the value ofNerwhich is directly proportional to the exact solutions is decreased with increasing the value ofk2.


## 6. Conclusion

Fuzzy theory provides a suitable way to objectively account for parameter uncertainty in models. Fuzzy logic approaches appear promising in preclinical applications and might be useful in drug discovery and design. Considerable progress has been made in the last few years in the development of computational approaches for prediction of drug absorption, distribution, metabolism, and excretion. Whilst several approaches have been developed in pharmacokinetics, most of these approaches have not yet been adequately used in the complex process such as prediction of metabolism, and they require further improvement.In summary in this research, a tau method based on the Jacobi operational matrix was utilized to numerically solve the PKPD equation, arising from drug assimilation into the bloodstream. The comparison of the results shows that the present method is a powerful mathematical tool for finding the numerical solutions of a generalized linear fuzzy fractional PKPD equation.Although we concentrated on applying our algorithm to solve fuzzy fractional PKPD equation, we show that such algorithm can be applied to solve other types of fractional equations models in science and engineering fields. Our algorithm for the fuzzy fractional PKPD equations is efficient and numerically stable. Numerical results are presented which exhibit the high accuracy of the proposed algorithm.By using the obtained operational matrix, the following objectives were achieved.(i)
The given fuzzy fractional differential equation was converted into a fuzzy algebraic system of equations, simplifying the solution procedure.(ii)
The method is computer oriented; thus, solving fuzzy fractional differential equations of different fractional ordersvbecomes a matter of changingvonly.



---

*Source: 304739-2013-12-15.xml*