*Department of Physics, Nigerian Defence Academy, P. M. B. 2109, Kaduna, Nigeria. ABSRTACT A first order stochastic semi-empirical model for pharmacokinetics is presented and the real response of drug concentration to vital pharmacokinetics parameters studied. By invoking Gaussian kinetics and the residual drug concentration eliminated, the probability densities and the response of concentration profiles are theoretically simulated, using empirical data based on our experience. The drug is administered for 3 days at regular time intervals of 3hr and 6hr, respectively, by refreshing the drug half-life.

Results show that the amount of drug residue decreases with increasing dose, but increases with increase in ingestion time interval for corresponding dose. It is also shown that the real drug concentration increases to a threshold and decreases marginally for subsequent dose. However it is difficult to predict the response of drug concentration with changes in ingestion time interval. We recommend that for higher drug concentration the half-life be increased. Our simulation results qualitatively agree with those documented in the literatures.

Keywords: Real drug concentration; Residue; first order pharmacokinetics; Gaussian kinetics; drug half-life; stochastic dynamics. 1. INTRODUCTION Pharmacokinetics is a branch of pharmacology dedicated to the study of time course of substances relationship with an organism. In practice it is applied mainly to drug substances. It concerns itself with all manner of compounds residing within the organism or system, such as nutrients, hormones and toxins, while pharmacodynamics explores what the drug does; pharmacokinetics considers what the body does with the drug.

Authorized drugs are those drugs commonly used by people for the treatment to their illnesses. The human body passes through various stages to process and distribute both authorised and prohibited drugs. These stages involve drug absorption into the body, drug distribution to various tissues and then the elimination or excretion (Pratt and Taylor, 1990; Amdur and Klaassen, 1991). Absorption is the first process when we take in drugs. This is the process in which the drugs pass through the administration into the bloodstream. Blood transmit the drugs to the various parts of the organ of the organ.

In this stage all the nutrients are absorbed by the blood and may have good or adverse effects to the body (Ellenhorn and Barceloux, 1998). The next stage is the distribution of the drugs. The blood is responsible for transporting the drugs to the organs in the body. The heart, the kidneys and the brain are organs in the human system that contains large capacity of blood. Due to this they are that main transporters of drugs in the human system. The last stage is the elimination. Drugs can be eliminated, as residue, through urinating, faeces, liver and even the mother’s milk.

These channels of drug excretion are essential in the execution of our body system. Drugs help in repairing the strength of an individual. Unfortunately, those harmful drugs may lead to various illnesses resulting into death. Prohibited drugs affect the brain and the body mechanisms (Pratt and Taylor, 1990). Patients commonly receive two or more drugs concurrently and most individuals who abuse drugs are poly-drug users. Multiple drug use may result in drug interactions. This occurs when the pharmacokinetics or pharmacodynamics of one drug is altered by another.

This concept is important to consider because interaction may result in decreased therapeutic efficacy or increased risk of toxicity. The degree of drug interaction depends on the relative concentrations and therefore dose and time (Roland and Tozer, 1989). Blood plasma is an essential fluid in understanding the transport of drugs in the human body. The packets transported by the human circulatory system are mainly contained in the blood plasma. This yellow liquid component of the blood constitutes about 55% of the total blood capacity.

It contains dissolved food substances and also serves as a medium for excretory product transportation. The blood flow is very important or decisive part of the human circulatory activities. Blood is transported to all the body tissues through an intrinsic network of blood vessels (Salloum, et al. , 2005). In the body tissue, the blood exciting the arteries and flowing into the capillaries is divided into blood in the core and blood flowing into other peripheral part of the human body such as the skin. Among the dissolved contents in the blood plasma, distribution of drugs is of utmost importance.

The mathematical treatment of drug distribution requires a thorough knowledge of the pharmacokinetics profile (Derendorf and Hochhaus, 1995). Along with the experimental approach, mathematical modelling has become a popular tool for analysis of drugs and cardiovascular systems (Eun et al. , 1995).

Global Dynamic Model (Gordon, 1959) has been developed for cardiovascular system; the Lumped Parameter Model (Hardy et al. , 1982; Masunzawa et al. , 1992) have been employed for controlling heart mechanisms and vessel hemodynamic; a model- based understanding of the fundamental of blood flow dynamics, or hemodynamic (Eun et al., 2004) is important for the diagnosis and management of diseases of the cardiovascular system, including coronary artery and heart muscle dysfunction, vascular disorders, and pulmonary disease; the One-and Two-Compartments Models, with advances in computer software, have been described (Hagan, 1996);

Physiological Models, in contrast to the Compartment Models, have been proposed (Rowland and Tozer, 1989); linear system theory has been employed to model many pharmacokinetics systems (Guyton et al. , 1984); the Capillary System Model was used to predict the blood flow in the arteries and veins (Salloum et al., 2005).

HPLC and fluorescence detection (Samanidou, et al. , 2005) have been employed in the simultaneous determination of quinine and Chloroquine by carrying out statistical evaluation with high accuracy.

The concept of zero or first order kinetics may be utilized to describe any rate process in pharmacokinetics (Pratt and Taylor, 1990). Therefore if we are discussing drug absorption, a drug exhibits zero order kinetics if a constant amount of drug is absorbed regardless of the dose. Conversely a drug exhibits first order absorption kinetics if the amount absorbed depends on the dose.

Most drugs exhibit first order elimination kinetics in which a constant fraction of drug is eliminated per unit time. This paper adopts the features of the one-compartment model (Hagan, 1996) and is somewhat in close comparison with the documented model (SDI Review team, unpublished).

Unlike the previous models we explore the dynamic behaviour of the drug concentration equation commonly derived in the literatures. The uniqueness in the present formulation is the elimination of the residue, using mathematical tools. This yields the real concentration of the drug.

In the frame of stochastic dynamics, we firstly formulate the probability density for the residue, which facilitates the derivation of the probability density for the real drug concentration as a complementary function. The isolated singularity is explored, which gives the residue, using the well known residue formula adopted from complex variable theory. We invoke Gaussian kinetics, based on the assumption of the Central Limit Theorem (Harry and Steven, 1994). Such an assumption allows the formulation of the first-order stochastic model from which the drug concentration profile evolves.

Our model permits theoretical simulation of the real response of drug concentration to vital pharmacokinetics parameters. The remaining part of this paper is structured as follows: in section 2 the probability densities and the stochastic model are formulated. The simulation results are discussed in section 3 and conclusion made in section 4. 2. MODEL FORMULATIONS.

2. 1THE PHYSICAL MODEL Similar to the One-Compartment model (Hagan, 1996) in which the entire human body is considered a single unit, the present model further describes the various stages of the drug dynamics in the body volume as shown in Fig.

1, in comparison with the One- Compartment model (Hagan et al. , 1996). We assume that the ingestion and the absorption of drug are contemporaneous and that the ingestion is at regular interval of time I hours, in doses of quantity C0, which is also the initial concentration, for a period of time. 2. 1. 1 THE ABSORPTION STAGE In this stage the initial concentration of the drug sample C0 is taken to be at time t=0.

The drug is absorbed at time t>0, the drug sample begins to react and its rate of reaction per time is the concentration C(t). 2. 1. 2 THE DISTRIBUTION STAGE.

We consider that the drug sample is distributed through the entire body as a single channel. At time t>0, the total concentration of the drug obeys the linear super-imposition: Ct=C0+CnI, (1) where n denote a particular dose count. After a certain time, the drug begins to decay. Salloum et al. (2005) considered that the drug is distributed at different channels in the system, contrary to the One- Compartment model. 2. 1. 3 THE ELIMINATION STAGE At this stage, due to certain processes (Derendorf and Hochhaus, 1995), which include elimination, the concentration is the sum of the pure and residual samples.

This residue is eliminated from the system. The concentration of the sample to the entire system to carry out certain action is thus a stochastic process. It is thus necessary to determine the probability density for the distribution of the residue. We try, in this paper, to derive the theoretical probability density which facilitates further mathematical derivations. Firstly we impose the mutual exclusiveness axiom: Pr+Pp=1 (2) Here Pr and Pp denote the probabilities of residue and dissolved drug samples respectively.

For preference we define the dissolved drug sample as the ‘real’ drug concentration, in our nomenclature. 2. 2. MATHEMATICAL FORMULATIONS 2. 2. 1ANALYSIS OF RESIDUE Consider a drug sample of known half-life T1/2 given every interval of time I in doses of quantity C0, for an extended period of time. For an initial dose C0 the concentration Ct of drug at any time t>0 is known to be described by the popularly known first order ordinary differential equation written as dC(t)dt=-kCt, t>0 (3a) C(0)= C0, (3b) where k is the elimation time rate or decay rate of drug sample.

The analytic solution of this equation can be straight forwardly be obtained by direct integration. This yields the solution; Ct=C0e-kt (4) At time t=I, the second dose of C0 is taken, increasing the concentration to CI=C0(1+e-kt) (5) The drug level immediately begins to decay. The mathematical expression due to this decay becomes dC(t)dt=-kCt, t>0 (6a) CI=C0(1+e-kt) (6b) Similarly the solution to this initial value problem (IVP) is Ct=C0(1+e-kt)e-k(t-I) (7) The next dose is taken at time t=2I. This yields the particular solution C2I=C0(1+e-kt)e-kI (8a).

The level increases to Ct=C0(1+e-kt+e-2kI) (8b) For successive doses, after (N+1) ingestion, the drug concentration is CNI=C0(1+e-kt+e-2kI+…+e-NkI). (9) This is the general solution popularly found in the literatures. Equ. (9) can be written in compact form as CNI=C0(1-e-(N+1)Ik)1-e-Ik (10) By partial fraction decomposition, this equation can be expressed as CNI=C01-e-Ik+C0e-(N+1)Ike-Ik-1 (11) In the limN>?

CNI we obtain the saturation concentration of drug given as Cs=C01-e-Ik (12) It is needful to remark that, albeit this saturation term has been given much attention in the literatures, the other term has never been explored.

Interestingly, the mathematical formula, Equ. (10), permits the synthesis of the residue. Denote this term as Cr. Thus Cr=C0e-(N+1)Ike-Ik-1 (13a) For successive dose we can express this term as a series of the first n terms Cr=n=1NC0 e-nIk(e-Ik-1) (13b) For the moment we conjecture that the residue is insitu in Equ. (13) and later justify that this residue is non-negative. It has been predicted (SDI Review team, suggested model) that for increasing time decay rate, larger than the injection time interval, it is possible that the residue of the drug increases.

More interesting is that analysis of drug in the human body needs to include the residue. In order to extract this residue, it is naturally assumed that the time course of drug (in the plasma) is linear.

We admit this deduction and employ the linear transformation; enx? 1+nx (14) This transformation reduces Equ. (13) into Cr=-n=1NC01+nkI(kI) (15) Here the accented variableCr denote the linearized form of Cr . In the literature (CRC Press, 1998) a plot of lnCNI against t gives a linear graph of the form Ct=? CnI+? , (16) Where the slope ? coincides with lnC0 and ? is the intercept.

In the realm of stochastic dynamics we define the moment generating function (MGF),? (t)l, in terms of ? and ?. Thus ; ? (t)l,=EetC=e? In? (? In) (17a) For generality we can define the characteristic function in terms of the set C as ? (t)C=E(e? ), (17b) where ? is a complex number given as; ?=inIC(nI) (17c) To preserve the linearity of the model, we freely admit that the characteristic function be defined by; ? (t)C=1, t? 0 (but zero else where) (18) This characteristic function is immediately invoked into Equ. (15) as follows: Cr=-n=1NC01+nkI(kI)? (t)C (19) To allow the convergence of Equ.

(19), a transformation of the form is necessary, as follows: Cr=-n=1NC0Ik2(n+1Ik)? (nI)C (20) Here the convergence is achieved such that 1Ik decreases monotonically. It is compelling to remark that the regularity in the injection time interval is a special case of the model described in the work (SDI Review team, suggested model). In their model, we appreciate that for equal time interval, the Heavy-side function invoked physically describes an integrated sharp impulse of unit magnitude. This is in coincidence with the form of characteristic function invoked in the present work.

As an implication, Equ. (20) shows that Cr has an isolated singularity at n=-1Ik with a pole of order 1. Following the well-known standard complex variable theory, we can calculate the residue,R, by the formula; R=limn>-1Ik{dp-1dnp-1n+1IkCrnI}P, (21) Where p is the order of the pole. Here p=1. We can immediately calculate the concentration of residue from this formula for the nth dose as; Cr(nI)=C0(Ik)2;t? 0 (22) 2. 2. 2 THE PROBABILITY DENSITY FUNCTIONS At any time t the theoretical probability function for the residue can be defined as Pr=C0Cs(nIk)2;t? 0 (23) Following the axiom, Equ.

(2), the probability function for the real drug concentration can be obtained as; Pp=1-C0Cs(nIk)2; t? 0 (24) We now invoke Gaussian kinetics thereby permitting the approximation of the probability of real drug concentration to Gaussian distribution, based on the Central Limit Theorem. This approximation is motivated by our conjecture that the eliminated residue is interpreted, in stochastic terms, as decreasing the skewness of the distribution. Recall that the normalized Gaussian distribution with flunctuation ? =1, can be written as; fCt=entire space12?

e-12(Ct)2=1 (25) In the realm of our earlier conjecture, the probability function, Equ. (24), can be normalized as; entire space(1-C0Cs(nIk)2)=1 (26) 2. 2. 3 THE STOCHASTIC MODEL Equating the last two equations and manipulating the resultant equation yields; Ct=ln(12? )-ln(1-C0Cstk2) (27) To formulate the first order stochastic model, Equ. (27) is simply differentiated once with respect to time and setting T1/2=ln2k (Hagan, 1996) we obtain; dC(t)dt=-2? IC0(T1/2)2(ln2t)2Cs-C0(T1/2)2t (28) This equation is the differential equation for the real concentration of drug, indexed in time. 3. RESULTS AND DISCUSSION.

3. 1RESULTS A numerical algorithm was employed to compute the concentration profiles, formulated in Equ. (28), using the familiar Euler’s method. The advantage of this method is due its simplicity in solving first order boundary value problems, in comparison with other existing numerical methods. The formula for the algorithm is given as: Ct+1=Ct+-2? IC0(T1/2)2(ln2t)2Cs-C0(T1/2)2t (29) The probability densities for the residue and the dissolved drug samples, as well as the numerical values of drug concentration profiles are calculated and the obtained results presented graphically (Fig.1-6).

The simulations are carried out theoretically by using the following data, based on our experience: C0=190mgkg, I=3hr, 6hr. The simulations are refreshed, for each value of I, using the following drug half lives: T1/2=7. 5hr, 15hr, 30hr, 60hr. All calculations are considered for 3 days of drug administration. In order to circumvent possible discontinuity in the numerical result att=0, we assumed thatC0=C1. The values of Cs, shown on table 1, were calculated using Equ. (12). Care must be taken that the saturation concentration does not imply the threshold concentration referred herein.

To clarify this, one may consider that while the saturation term refers to the global optimum of the drug concentration the threshold term refers to the local optimum of the drug concentration. It is worth remarking that the values of the first dose probabilities were taken to be zero, to relegate spurious values. Albeit, this has no effect on the overall simulation result. The simulations for I=3hr are labelled A1-A4 while for I=6hr are labelled B1-B4 for clarity (Table 1). Table 1: Summary of simulation data I(hr) label T12(hr) C0(mgkg) Cs(mgkg) Duration (days) A1 7. 5 190 783. 58 3 3 A2 15 190 1468. 05 3 A3 30 190 380.

06 3 A4 60 190 5571. 05 3 B1 7. 5 190 446. 43 3 6 B2 15 190 784. 48 3 B3 30 190 1468. 32 3 B4 60 190 6052. 41 3 4. 2. DISCUSSION The reliability of our model can be seen already in the simulations of the probability densities (Tables 2 and 3) for residues, the probability densities for real drug concentration (Fig. 2 and 3), and more importantly the real drug concentration profiles (Fig. 4 and 5). Tables 2 and 3 show that the residues decrease with increase in dose , but increases with increase in drug half-life. Comparing the figures, it can be observed that the residue is higher for larger ingestion time interval.

This may be likened to susceptibility of drug to other foreign substances due to prolonged decay time. This is in agreement with the influence of active metabolites interacting with the parent drug (CRC Press, 1998; Rowland and Tozer, 1989). It is simply expected that the corresponding probabilities (Fig 2 and 3) for the real drug concentration show complementary features to their respective residues. However one can probably deduce that the amount of residue becomes insignificant for higher dose at regular interval.

The real drug concentration profiles (Fig.4 and 5) show significant results. The drug concentration shows initial decline after peaking. The threshold time (time-to-peak) can be seen to vary with drug half-life and ingestion time interval. These are in good agreement with earlier findings (e. g. , Hagan 1996; Rowland and Tozer 1989). Rowland Tozer (1989) confirmed the decline in drug concentration as drug tolerance which is influenced by ingestion time. Our simulations also confirm that the degree of tolerance is not consistent with drug half-life. This can be easily observed by the inconsistency shown (Fig. 4 and 5) for T1/2=30 hrs.

This inconsistency infers that drug concentration is sensitive to the correlation between drug half-life and ingestion time-interval. Our simulations are not consistent with response real drug concentration to variation in ingestion time interval (Fig. 4 and 5). This can be observed by comparing the numerical values for I=3hr and I=6hr (Fig 4 and 5) at corresponding ingestion time. Fortunately it has been documented in the literature (CRC Press, 1998) that such difficulty may be likened to factors such as irregular change in blood pressure with time which possibly results in misinterpretation of clinic picture.

However it is probable that the difference in drug concentration at different ingestion time is insignificant for very large half-lives. Our numerical findings are less consistent with models simulating response of drug concentration at unequal ingestion time interval (e. g. , SDI Review suggested model). It would be interesting to understand what assumptions of these models differ. Finally our model is amenable to simulate response of free drug concentration of specific drugs with known empirical data. Table 2. Probability values for residue sample at 3hr ingestion interval N| Ingestion time (hr)| Probability values at half life 7.

5hr| Probability values at half life 15hr| Probability values at half life 30hr| 1| 3| 0. 000| 0. 000| 0. 000| 2| 6| 0. 789| 0. 000| 0. 000| 3| 9| 0. 351| 0. 748| 0. 000| 4| 12| 0. 197| 0. 421| 0. 871| 5| 15| 0. 126| 0. 269| 0. 557| 6| 18| 0. 088| 0. 187| 0. 387| 7| 21| 0. 064| 0. 137| 0. 248| 8| 24| 0. 049| 0. 105| 0. 218| 9| 27| 0. 039| 0. 083| 0. 172| 10| 30| 0. 032| 0. 067| 0. 139| 11| 33| 0. 026| 0. 056| 0. 115| 12| 36| 0. 022| 0. 047| 0. 097| 13| 39| 0. 019| 0. 040| 0. 083| 14| 42| 0. 016| 0. 034| 0. 071| 15| 45| 0. 014| 0. 030| 0. 062| 16| 48| 0. 012| 0. 026| 0. 054| 17| 51| 0. 011| 0. 023| 0.

048| 18| 54| 0. 010| 0. 021| 0. 043| 19| 57| 0. 009| 0. 019| 0. 039| 20| 60| 0. 008| 0. 017| 0. 035| 21| 63| 0. 007| 0. 015| 0. 032| 22| 66| 0. 007| 0. 014| 0. 029| 23| 69| 0. 006| 0. 013| 0. 026| 24| 72| 0. 006| 0. 012| 0. 024| | | | | | Table 3. Probability values for residue sample at 6hr ingestion interval N| Ingestion time (hr)| Probability values at half life 7. 5hr| Probability values at half life 15hr| Probability values at half life 30hr| 1| 6| 0. 000| 0. 000| 0. 000| 2| 12| 0. 346| 0. 785| 0. 000| 3| 18| 0. 154| 0. 350| 0. 749| 4| 24| 0. 087| 0. 197| 0. 421| 5| 30| 0. 055| 0. 126| 0.

270| 6| 36| 0. 039| 0. 088| 0. 187| 7| 42| 0. 028| 0. 064| 0. 135| 8| 48| 0. 022| 0. 049| 0. 105| 9| 54| 0. 017| 0. 039| 0. 083| 10 | 60| 0. 014| 0. 032| 0. 067| 11| 66| 0. 011| 0. 026| 0. 056| 12| 72| 0. 010| 0. 022| 0. 047| 13| 78| 0. 008| 0. 019| 0. 040| Figure 2: Probability densities for real drug concentration profiles against time for ingestion time interval of 3hr Figure 3: Probability densities for real drug concentration profiles against time for ingestion time interval of 6hr Figure 4: Real drug concentration (mg/kg) profiles against ingestion time for ingestion time interval of 3hr.

Figure 5: Real concentration of drug (mg/kg) against ingestion time (hr) for ingestion time interval of 6hr 4. CONCLUSIONS Validity of our model for first order decay pharmacokinetics has been shown. The probability densities and the real response of drug concentration to drug half life and ingestion time interval have been theoretically studied. Irrespective of the half life, all simulations have given consistent and reliable results. The eliminated residue has significant influence on the drug concentration profile, but the quantitative difference in the values of the drug concentration is not very important.

The major discrepancy lays in the sensitivity of drug concentration profiles to half life. It is thus necessary to analyse the drug residue. This gives quantitative information about the real concentration of the drug in the human body. We admit that the use of such mathematical models should not be treated in absolute terms. They are only attempts to understand the behaviour of the drugs in the human body. ACKNOWLEDGEMENT The authors wish to acknowledge the SDI review team for suggested model and constructive criticism. REFERENCES.