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

and Applied Analysis 3 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 PKPDmodels. 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 Section 3, 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. Let u be a fuzzy set in R. u is called a fuzzy number if: (i) u is normal: there exists x 0 ∈ R such that u(x 0 ) = 1; (ii) u is convex: for all x, y ∈ R and 0 ≤ λ ≤ 1, it holds that u (λx + (1 − λ) y ≥ min {u (x) , u (y)}) ; (1) (iii) u is upper semicontinuous: for any x 0 ∈ R, it holds that u (x 0 ) ≥ lim x→x ± 0 u (x) ; (2) (iv) [u]0 = supp(u) is a compact subset of R. In this paper, the set of all fuzzy numbers is denoted by RF. Definition 2. Let u ∈ RF and r ∈ [0, 1]. The r-cut of u is the crisp set [u]r that contains all elements withmembership degree in u greater than or equal to r; that is [u] r = {x ∈ R | u (x) ≥ r} . (3) For a fuzzy number u, its r-cuts are closed intervals in R, and we denote them by [u] r = [u r 1 , u r 2 ] . (4) According to Zadeh’s extension principle, the operation of addition on RF is defined as follows: (u + V) (x) = sup y∈R min {u (y) , V (x − y)} , x ∈ R, (5) and scalar multiplication of a fuzzy number is given by (k ⊙ u) (x) = { { { u( x k ) , k > 0,


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 "pHpartition" hypothesis to predict the fraction of dose absorbed [3][4][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][9][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].

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][26][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][30][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][39][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][44][45][46], Chebyshev polynomials [47][48][49], Laguerre polynomials [50][51][52][53], and Bernstein polynomials [54][55][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][62][63][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][65][66][67][68][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 Section 3, 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.

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.

The Fuzzy Settings Definitions
Definition 1.Let  be a fuzzy set in R.  is called a fuzzy number if: (i)  is normal: there exists  0 ∈ R such that ( 0 ) = 1; (ii)  is convex: for all ,  ∈ R and 0 ≤  ≤ 1, it holds that (iii)  is upper semicontinuous: for any  0 ∈ R, it holds that In this paper, the set of all fuzzy numbers is denoted by R F .Definition 2. Let  ∈ R F and  ∈ [0, 1].The -cut of  is the crisp set []  that contains all elements with membership degree in  greater than or equal to ; that is For a fuzzy number , its -cuts are closed intervals in R, and we denote them by According to Zadeh's extension principle, the operation of addition on R F is defined as follows: and scalar multiplication of a fuzzy number is given by where 0 ∈ R F .
Definition 3 (see [70]).The distance (, V) between two fuzzy numbers  and V is defined as where is the Hausdorff distance between []  and [V]  .
It is easy to see that  is a metric in R F and has the following properties (see [70,72]): Definition 4 (see [73]).Let  and  be the two fuzzynumber-valued functions on the interval [, ], that is, ,  : [, ] → R F .The uniform distance between fuzzy-numbervalued functions is defined by Remark 5 (see [73]).Let  : [, ] → R F be fuzzy continuous.Then from property (iv) of Hausdorff distance, we can define Definition 6 (see [74]).Let ,  ∈ R F .If there exists  ∈ R F such that  =  ⊕ , then  is called the H-difference of  and , and it is denoted by  ⊖ .
In this paper, the sign "⊖" always stands for H-difference, and note that  ⊕  ̸ =  + (−).Also throughout the paper, it is assumed that the Hukuhara difference and generalized Hukuhara differentiability exist.
The matrix form of the above equations is where the coefficient matrix  = (  ), 1 ≤ ,  ≤  is a crisp  ×  matrix and   ∈ R F , 1 ≤  ≤ .This system is called a fuzzy linear system (FLS).
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].
Let  be the smallest integer that exceeds V, then Caputo's fractional derivative operator of order V > 0 is defined as where For the Caputo derivative, we have Caputo's fractional differentiation is a linear operation, namely, where  and  are constants.
Theorem 10 (generalized Taylor formula [78]).Suppose that    () ∈ (, ] for  = 0, 1, . . .,  + 1 where 0 <  ≤ 1; then with where Here,    shows fractional derivative operator in the Caputo sense, and = is called fractional Taylor expansion with degree  of , and    (, ) is reminder term in fractional Taylor expansion of .Now, we present some fuzzy fractional notations which are used later in the paper.
indicates the set of all fuzzy-valued functions which are continuous up to order .
(iv)  R F [, ] 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.
The Riemann-Liouville integral of fuzzy-valued function  is defined as Definition 12 (see [63]).
We say that () is fuzzy Caputo fractional differentiable of order 0 < V ≤ 1 at  0 , if there exists an element or or or For the sake of simplicity, we say that the fuzzy-valued function  is  [(1) − V]-differentiable if it is differentiable as in Definition 13 case (i), and  is  [(2) − V]-differentiable if it is differentiable as in Definition 13 case (ii), and so on for the other cases.

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.
Drug intake GI tract Digestion Blood Tissue For each compartment by applying the balance law we can obtain 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.

Case of a Single Cold Pill.
Let us consider () to be the amount of drug in the GI-tract at time  and (), the amount in the bloodstream at time .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], where  0 is the amount of a drug in the pill and  1 is a positive coefficient.We suppose that (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, with  2 another positive constant.Decongestant and an antihistamine are the component of the cold pill, and the coefficient of proportionality ( 1 and  2 ) 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 factor   1  , integrate, and then use the initial data to obtain Insert this formula for () into the second rate equation, which becomes In this paper, we are concerned with fractional time derivatives regarding (17) and (18).Having defined   V , 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, ], 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 point  = 0 [81], which are functions of .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: In the standard direct action model, the effect at time , (), is expressed by an arbitrary (memory-less) function of drug concentration in the effect site at time , (()); however to generate a wider class of relationships, we assume that the effect at time  ∈ [0, 1] is related to the fuzzy Caputo fractional derivative of ().So, we have ] is a continuous fuzzyvalued function and   V 0 + denotes the fuzzy Caputo fractional derivative of order V ∈ [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][58][59].
We have assumed that  1 ̸ =  2 , 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 are  1 = 0.6931 (hour) −1 and  2 = 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 define  2 as clearance coefficient of medication from the blood.The value of  2 often 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 that  2 is much lower, which is normally observed in old and sick cases.
We investigate sensitivity of the medication over a 24hour period by keeping the value of  1 fixed at some value such as 1.386 but setting  2 = 0.01386, 0.06386, 0.1386, 0.6386, 1.386 ( = 1).Different values of  2 correspond 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 high  2 , 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 that  2 is constant ( 2 = 0.0231) and the value of  1 varies ( 1 = 0.06931, 0.11, 0.3, 0.6931, 1.0, and 1.5).In this case, the  1 is a constant for the GI tract which is analogous to the  2 for the bloodstream.

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 of  1 and  2 .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] where the fuzzy coefficients   are gained by and (,) 1 () = a  ⊙ Φ  () , (52) where the fuzzy shifted Jacobi coefficient vector   and shifted Jacobi polynomials vector Φ +1 () are presented by We can state the fuzzy approximate function (52) in the parametric form as follows.
Let us define the following notations: 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 Denoting we investigate that ( 58) is equivalent to the matrix system The elements of the matrices mentioned above are determined in the next theorem.
Theorem 21 (see [59]).Let us denote   = ( (V) 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.
Proof.Since  V 0+ () ∈ (0, ],  = 0, 1, . . .,  + 1,  can be expanded to the fractional Taylor series and its reminder term is Therefore, From ( 65) and triangle inequality, we obtain 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.
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.

Numerical Results
In this section, the fuzzy fractional PKPD model ( 49) is solved for different values of  1 and  2 by 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 [  ]  = [   ,    ] and the corresponding exact solutions (; ) = [(; ), (; )], that is, are considered.Now, we recall the FFDE (49) as follows: ] is a continuous fuzzyvalued function and   V 0 + denotes the fuzzy Caputo fractional derivative of order V ∈ [0, 1].
Two cases are considered for solving by the proposed technique.In the first case, we assume that  2 is unchangeable and different values of  1 are substituted in (49) to get the fuzzy approximate solution.Conversely, in the second case we try to solve the problem such that  1 has an invariable value and  2 is varied.The problem is analyzed for both cases in details.
Case I. Let us consider from Section 3 that  2 = 0.0231, but let  1 vary (e.g., 0.6931, 0.11, and 0.3).So (49), with the assumption, namely,  2 = 0.0231 and Dose  = 1, is as follows: By using Theorem 7 for the above equation under  [(1)−V] differentiability, we have the following systems: Solving ( 72) leads to determining the exact solution of (71) as follows: We seek the fuzzy approximate solutions by applying the technique described in Section 4 with  = 2 and as Here, we have  0.85 = ( 0 0 0 1.8639 0.3901 −0.1755 −0.3901 4.5267 0.8696 where  =  = 0. Regarding (54), we have By solving the fuzzy linear system (60) in the parametric form, one can get fuzzy unknown coefficients in (76) as where we assume that -cut = 0.1 in (77).The comparison between absolute errors of different  1 obtained 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 of  1 with the fractional order V = 0.85 in Figure 4.
The absolute error value for -cut varied from 0 to 1 for different value of  1 which is calculated in Table 1.As it can be observed at a constant -cut by increasing the value of  1 which 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 the  1 [12].It should be noted here that the value of    is directly proportional to the variation of exact solutions and corresponding fuzzy approximate solution.As a result, increase of    with increasing  1 is also expected.
Case II.In this case, we assume that  1 = 1.386 and is not variable, but  2 takes different values (i.e.,  2 = 0.01386, 0.06386, and 0.1386) and Dose  = 1.So (70) alters with, for example,  2 = 0.01386 as follows: with the exact solution as: By applying the JT method in Section 4, we can get the fuzzy unknown coefficients {  }  =0 by using (58) which is in the parametric form as follows: where  (V) is given as with  =  = 0.5 and  = 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.We compared JT method results for different values of  2 and 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  = 0 and  = 0. Additionally, we compared the approximate solutions obtained by the present method at  = 4, 7, 9, and 11 which 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 of  2 at V = 0.75.It can be seen that for all the values of  2 the approximate solution is the fuzzy number.
In Table 3, the value of absolute error for -cut varied from 0 to 1 for different value of  2 is calculated.At a constant -cut value, a descending trend is observed with increasing the value of  2 .This trend is opposed to the trend observed in Table 1 for different value of  1 . 2 is clearance constant and it is a coefficient for () which is the amount of drug in the bloodstream in (44).We already explained that by decreasing the value of  2 ,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 of  2 appears with negative sign in the main differential equation as we expected, the value of    which is directly proportional to the exact solutions is decreased with increasing the value of  2 .

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 orders V becomes a matter of changing V only.

Figure 1 :
Figure 1: Schematic of input/output compartment for drug assimilation.

Table 1 :
The results of the proposed method for Case I with V = 0.85,  =  = 0, and  = 8.

Table 3 :
The results of the proposed method for Case II with V = 0.95,  =  = 0.5, and  = 7.