A physiologically-oriented mathematical model for the description of in vivo drug release and absorption

This paper focuses on a physiologically-oriented mathematical model aimed at studying the in vivo drug release, absorption, distribution, metabolism and elimination (ADME). To this purpose, the model accounts for drug delivery from an ensemble of non-eroding poly-disperse polymeric particles and the subsequent ADME processes. The model outcomes are studied with reference to three widely used drugs: theophylline, temazepam and nimesulide. One of the most important results of this study is the quantitative evaluation of the interplay between the release kinetics and the subsequent ADME processes. Indeed, it is usually assumed that in vivo drug release coincides with in vitro so that the effect exerted by the ADME processes is neglected. In addition, the proposed model may be an important tool for the design of delivery systems since, through proper changes, it could also account for different oral delivery systems.


Introduction
Quite recently, an initiative of the Pharmaceutical Research and Manufacturers of America (PhRMA) tried to evaluate the reliability of predictive models in terms of of drug efficacy, safety and properties estimation [1][2][3]. For this purpose, 108 clinical compounds (22 % acids, 46 % bases, 18 % zwitterions and 14 % neutrals), supplied by 12 PhRMA member companies, have been considered. The overall goal of this study was to assess the predictability of human pharmacokinetics (PK) from preclinical data and to provide comparisons of available prediction methods from literature. In particular, Physiologically Based PK (PBPK) models, despite their limited application in drug discovery and development, due to their mathematical complexity and requirement of many experimental data, have been considered. Indeed, these models have doi: 10.5599/admet.2.2.34 81 become more attractive for the prediction of human PK profiles thanks to the prediction of hepatic metabolic clearance and tissue distribution from in vitro and in silico data as well as the advent of commercially available and user-friendly software [4]. PBPK models have showed, in 69% of the cases, a medium to high degree of accuracy in simulating the time course of the drug plasma concentration in human following intravenous administration. On the contrary, the accuracy decreases to 23% in case of oral administration [5]. Accordingly, the authors of this study concluded that i) PBPK models can be a useful basis for a more rational selection of first-in-human phase dose ranges and that ii) there is a need for better predictions of human pharmacokinetics following oral administration. Indeed, predicting the time course of drug concentration in plasma/blood following oral administration presents the not negligible difficulty of simulating both the delivery and the absorption stages which is obviously not needed with intravenous administration. When efforts have been focused on taking into account one of these two stages, they were directed to the absorption [6], neglecting the delivery stage. As a matter of fact, both the PhRMA study and common commercial software consider the release stage either by assuming a fixed release kinetics or by providing the model with an in vitro release profile assumed to be representative of in vivo release [5]. In this light and considering that oral administration is by far the most common route of drug administration [5], we believe it is worth working on the improvement of PBPK models devoted to simulate PK following oral administration. In addition, we are convinced that the new challenge in the drug delivery field, for mathematical modelling, is to combine of mechanistic theories able to realistically describe the simultaneous processes of drug release and the subsequent ADME processes within the human body [7,8]. Indeed, this step is a fundamental prerequisite to get the final objective of building up a mathematical model able to predict pharmacodynamics, i.e. the clinical effect of a drug released from a specific delivery system [9]. Nowadays, many powerful commercial software exist for the simulation of PK phenomena. Among them, Simcyp, Gastroplus and PK-SIM can be cited. These commercial PBPK models combine a variety of functions for defining release from oral formulations with sophisticated models for the description of ADME processes. In particular, Simcyp is able to account for population variability of parameters and their covariation where known. Our mathematical model, on the contrary, tries to establish a deep connection between the ADME phenomena, described by a recently developed simple PBPK model [10], and the drug delivery step, described by a model [11,12] that proved to be reliable. In particular, the considered delivery system consists of an ensemble of non-eroding poly-dispersed polymeric particles hosting the drug inside polymeric meshes. Particulate systems, indeed, present remarkable advantages over the single unit devices including easier dispersion within the stomach that can result in an appreciable reduction of local drug concentration, and reduction of the insurgence of gastric irritation [13]. Moreover, particulate systems are very versatile, one more reason for their wide use.

PBPK model
The mathematical representation of the human body is essentially based on the PBPK model presented by Di Muria and co-workers [10]. The gastrointestinal (GI) tract is described by a continuous hollow cylinder subdivided into three zones: the gastric lumen (GL) the small intestine lumen (SIL) and the large intestine lumen (LIL) (see Figure 1). The different physiological and morphological conditions (for example pH and internal radius) in the GI tract, account for different permeability (P GL , P SIL and P LIL , respectively), drug solubility (C sGL , C sSIL and C sLIL , respectively) and internal radius ( The GI tract is crossed by the release environment, an aqueous volume (V r ), containing the delivery system 82 (an ensemble of drug loaded, non-eroding, poly-dispersed polymeric particles), that, due to peristaltic movements, flows downwards from the GL to the SIL and eventually to the LIL. For modelling purposes, it is assumed that the release environment shape is that of a hollow cylinder whose external radius coincides depending on the GI position reached by the release environment. When, on the contrary, R i is > 0, the release environment assumes the shape of a hollow cylinder. In both cases, the length of the release environment is defined by (see Figure 1): Due to the well-known differences in the morphology/physiology of each zone, it is supposed that the release environment speed (v re ) can vary along the GI tract and this is indirectly considered for by setting three different residence (transit) times (t GL , t SIL and t LIL , respectively), one for each zone. Release environment velocity (v re ) and position in the GI tract (P os ) are connected by the following relation: P os = ∫ 0 t v re *dt, where t is time. Thus, assuming a constant v re , the transit time (t t ) of the whole GI tract is given by t t = P os /v re where, now, P os indicates the GI tract length. Alternatively, always assuming a constant v re , after time t, the position gained by the release environment is given by: P os = v re *t.
Once released, the drug can undergo elimination in the release environment (elimination constant RE el k , dimensionally time -1 ) and it can cross, due to passive diffusion, the GI mucosa according to the local permeability that depends on the V r position in the GI tract. Then, the drug reaches the gastrointestinal circulatory system (GICS), i.e. the ensemble of blood vessels around the GI tract. The compartment named GICS (Gastro Intestinal Circulatory System) is new to pharmacokinetic modeling, firstly introduced by [10]. Basically, it corresponds to the ensemble of the mesenteric artery, the portal vein and the ensemble of microcirculatory gastrointestinal vessels. This section of the whole circulatory system behaves differently from the remaining. In the case of parenteral administration, the mesenteric artery transports the drug towards the gastrointestinal tract. In the case of enteral administration, on the contrary, the portal vein transports the drug from the gastrointestinal tract towards the liver (where the drug experience the socalled first-pass effect) to finally get the remaining part of the circulatory system (here called "plasma") and the other tissues. The introduction of the GICS has allowed to effectively transform the very detailed pharmacokinetic model proposed by [14], (consisting in 21 compartments and 38 ordinary differential equations with roughly one hundred parameters), into a simple and effective model (7 compartments and 7 ordinary differential equations, with about twenty parameters). This effort saw several researchers spending lot of work, with limited and ineffective results (see, for example (see, for example, [15]). The GICS is assumed to be a well-stirred environment characterized by its own constant volume (V GICS ) and time dependent drug concentration (C GICS ).
The portal vein ensures the convective drug transport (blood volumetric flow Q PV ) from the GICS to the liver (volume V L ). In this well-stirred environment, drug concentration (C L ) varies due to hepatic clearance (C LH ; dimensionally a volumetric flow) and convective transport by the hepatic vein and the hepatic artery. The hepatic vein connects the liver to the highly perfused tissues and organs (where drug concentration is assumed equal to plasma's) by a volumetric blood flow Q HV , while the hepatic artery brings back the drug from the plasma environment to the liver (volumetric blood flow Q HA ). Once in the plasma, the drug can undergo elimination ( k must be volumetric flows [12]. However, as this is not a common choice and it can generate confusion in the reader, model equations have been re-arranged so that more common drug elimination (k el , time -1 ) and plasma -tissues exchange constants (k PT , time  (8)) that expresses the ratio between the plasma (V P ) and the tissues (V T ) volumes. Accordingly, the PBPK model discussed and reported in Figure 1, is mathematically represented by the following system of ordinary differential equations where t is the independent variable, time: where: Eq. (2) allows the computing of the drug concentration in the GICS, Eq. (3) evaluates the liver drug concentration, Eq. (4) determines the evaluation of the plasma drug concentration, Eq. (5) calculates the tissues drug concentration and, finally, Eq. (6) allows calculation of the amount of drug eliminated in the release environment, in the liver, and in the plasma. Inspection of Eq. (2) reveals that the diffusive drug flux from the release environment to the GICS requires knowledge of the drug concentration in the release environment (C r ) that is evaluated in the following section.

Delivery model
Briefly, our model assumes that the drug release kinetics are essentially determined by (i) polymeric particles size distribution in the dry state (ii) particle swelling upon contact with the release environment fluid (solvent) (particle erosion is considered negligible, this being true for highly crosslinked polymeric networks) (iii) drug dissolution and diffusion inside particles and (iv) possible drug re-crystallisation upon contact with the release environment fluid. The first point is achieved by assuming that dry particle size distribution can be conveniently represented by a Weibull distribution function [16]: where R p and R min are, respectively, the generic particle radius and the minimum particle radius,  and  are two parameters regulating the Weibull size distribution, while V 0 and V arethe total volume occupied by the ensemble of polymeric particles and the volume occupied by particles having a radius lower than or equal to R p , respectively. In order to account for particles polydispersity, the Weibull size distribution (Eq. (9)) is subdivided into N c classes (from i = 1 to N c ) on the basis of the particles radius. Accordingly, particles of class i th share the same radius R pi . For each class, the swelling process can be evaluated by the following mass balance referred to the swelling agent (solvent): where J fi and J ri indicate the Fickian and non-fickian component of the solvent flux, respectively. J ri represents a delayed flux induced by the polymer/solvent viscoelastic character. D 0 and D r are, the Fickian and the non-Fickian solvent diffusion coefficient, respectively, while  is the polymer-solvent relaxation time: where C peq is the solvent concentration in the completely swollen polymeric network (thermodynamic equilibrium),  eq and D eq are, respectively, the relaxation time and the solvent diffusion coefficient in equilibrium conditions while g and f are two model parameters evaluated in [11] by data fitting. The existence of two contributions in the solvent flux is motivated by the fact that the "polymer + solvent" system is viscoelastic, i.e. polymeric chains undergo rearrangements that develop with time (relaxation) upon solvent uptake. In order to account for this phenomenon, Camera-Roda and Sarti successfully proposed the approach expressed by Eq. (11) [17]. Eq. (10) can be solved assuming that, initially, particles are solvent-free, that solvent flux, for symmetry reasons, is zero in the particle centre and that solvent concentration at the particle surface (C pint ) is defined by: The swelling part of this model is completed by assuming that the solvent-polymer mixing is ideal, i.e., upon swelling, particle volume can be simply computed anytime by summing polymer and solvent volume. Inspection of Eqs. (10) -(13) clearly reveals that solvent diffusion is assumed independent of the presence of drug, this being rigorously true only for low drug concentrations in the particles.
Drug transport in the polymeric network (point iii) is modelled according to the following equations: Eq. (14) represents the drug mass balance relative to the "i th " class of particles, i.e. the class of particles characterised by radius R pi . This equation implies the existence of dissolved and undissolved drug inside the swelling polymeric network. While undissolved drug cannot move, dissolved drug mobility is controlled by its diffusion coefficient D that depends on the local solvent concentration C pi according to the Peppas -Reinhart theory [18] (see Eq. (15)). In Eq. (15), D ds represents the drug diffusion coefficient in the pure solvent while  i is the local polymer volume fraction and  s indicates solvent density. While C i is the position-dependent concentration of the dissolved drug, am di C , nc di C and mc di C are, respectively, the positiondependent concentrations of the undissolved drug in the amorphous, nano-crystalline and microcrystalline state. Indeed, in polymeric matrices, depending on the drug loading technique adopted, the drug can be found in amorphous, nanocrystalline and microcrystalline state [16]. As nano-crystals and amorphous drug are characterised by a higher solubility with respect to that of microcrystals [19], their presence can have a significant effect upon drug bioavailability. Whatever the state is, the drug dissolution rate (Eqs. (16), (18), (19)) is assumed to be proportional (K being the proportionality constant) to the difference between the local drug concentration C i and drug solubility ( am s C , nc s C , mc s C for amorphous, nanocrystal and microcrystal drug, respectively). In principle, the model could also account for the pH dependence of am s C , nc s C , mc s C . This, of course, would require the definition of a proper mathematical function linking solubility and pH. As the amorphous drug is often not stable and, upon dissolution, generally transforms into the more thermodynamically stable microcrystal condition (re-crystallisation; point iv)), its solubility reduces with time according to Nogami theory [20] (see Eq. (17)) where K r is the recrystallisation constant and am s0 C is the maximum amorphous drug solubility. Obviously, drug recrystallization can also take place in the release environment and Eq. (20) serves to account for this phenomenon. In particular, M c and K rb are, respectively, the amount of re-crystallized drug up to time t and the re-crystallization constant in the release environment. The solution of Eq. (14) is performed by assuming that, initially, no dissolved drug exists inside the particles (only solid or amorphous drug can be found in the polymeric network) and that the solid and/or amorphous drug concentration are uniform inside the particles (the model could easily account for more complex drug distributions). Finally, for doi: 10.5599/admet.2.2.34 87 symmetry reasons, drug flux in the particle center is set to zero and the usual drug partitioning condition at the particle/solvent interface is considered (K p is the drug partition coefficient).
In order to deal with the balance between problem unknowns (C GICS , C L , C p , C T , C Pi , C i , C r ) and equations, it is necessary to add a further equation for the definition of C r . In particular, the overall drug mass balance made on drug is considered: where M 0 is the administered drug dose, N c is the number of classes into which the continuous particle size distribution has been subdivided into (Eq. (9)) and N pi is the number of particles belonging to the "i th " class (N pi can be simply evaluated on the basis of Eq. (9)). Eq. (21)  Due to the simultaneous presence of ordinary and partial differential equations, the solution of the model calls for a numerical algorithm. In particular, the ordinary differential equations have been converted into linear algebraic equations according to the implicit Euler method [21] while the implicit control volume [22] (10 particles classes, 15 radial subdivisions) has been considered for the partial differential equations conversion into linear and non linear algebraic equations. In order to ensure numerical solution accuracy, the integration time step was set to 1 s. The relaxation method [21] has been used for the iterative solution of our algebraic equations system. The solution of each iteration was carried out according to the LU decomposition approach [21]. A home-made FORTRAN code was built up for the numerical solution.

Results and Discussion
One of the critical factors connected to the use of PBPK models is the availability of reliable values for the model parameters. Of course, this is not a trivial problem [1][2][3][4][5] because different values for the same parameter can be found in the literature due, in some cases to considerable inter-subjects variability [23 -29]. Accordingly, in order to simplify comparison of the three different drugs considered in this paper (theophylline, temazepam and nimesulide), we decided to rely exclusively on the information reported in [10] as source of PK parameters values.
The first drug considered, theophylline, also known as 1,3-dimethylxanthine, is a methylxanthine drug that finds a wide clinical use in the treatment of respiratory diseases such as asthma. Its chemical/physical and PK characteristics are reported in Table 1. Figure 2A shows the results of model simulations when assuming an orally administered dose of 100 mg dispersed in a polymeric carrier with a volume of 1 cm 3 . In particular, the continuous thick solid line refers to theophylline concentration in the plasma (C p ) in the case of only one spherical matrix of 6200 m radius. The initial part of the curve, characterised by zero drug concentration, is related to drug residence in the stomach where permeability has been set equal to zero (P GL = 0). Then, C p slowly increases to a maximum followed by a subsequent drop off till the complete drug disappearance from the plasma.  The C p increase/decrease trend corresponds to the drug absorption in the small and large intestine and is the result of both the simultaneous action of drug release from the matrix and the following ADME phenomena. The main reason for the low C p values is the very slow drug delivery kinetics (low ratio between the release surface and the matrix volume, typical of big spheres) from the matrix as depicted in Figure 2B (see the continuous thick solid line; C r indicates the time evolution of drug concentration in release environment). By keeping constant the delivery system volume (1 cm 3 ), its splitting into all equal smaller particles leads to a much faster C p increase. In particular, the thin solid line and the dashed line of Figure 2A represent, respectively, the C p trend referred to particles characterised by a radius equal to 620 m (about 10 3 particles) and 62 m (about 10 6 particles). It can be seen that, in both cases, the drugconcentration rise is fast and C p reaches about 5 g/cm 3 after 3.5 h in the first case (620 m) and 6.5 g/cm 3 after 2.5 h in the second case (62 m). Further reductions of particles size do not reflect in an appreciable variation of the C p trend. Indeed, as shown in Figure 2B, when particles radius is around 62 m, the drug, totally, and very rapidly, dissolves in the release environment (stomach fluids) to get the maximum value compatible with the administered dose, i.e. 400 g/cm 3 . In this case, no control on drug release is exerted and the delivery system simply ensures a prompt and total release in the stomach. In other words, no differences occur in the C p trend induced by this formulation or by a drug solution of the same dose. It is worth mentioning that the C p trend of the 62 m particles can be approximately obtained by adopting an ensemble of poly-dispersed particles described by eq.(9) considering the following parameters:  = 1000 m,  = 2, R min = 0 m and R max = 1200 m (see the grey line in Figure 2A). The difference of using a poly-dispersed distribution relies on a smaller drug concentration (C r , see grey line in Figure 2B) in the stomach fluids, which is advantageous in the case of stomach wall aggressive drugs.  89 Finally, it is worthwhile underlying that, based on the PK parameters set in Table 1, the time evolution of drug concentration in the CIGS (C GICS ) and in the liver (C L ) is substantially equal to that of C p . This depends on the small CIGS (V GICS = 10 cm 3 ) and liver (V L = 1500 cm 3 ) volume in comparison to the plasma (V P = 11000 cm 3 ) and tissues (V T = 24000 cm 3 ) volumes and the relatively high volumetric flow connecting the various compartments (Q HV , Q PV and Q HA ; see Figure 1 and Table 1). On the contrary, drug concentration (C T ) in the scarcely perfused tissues and organs (tissues in Figure 1) monotonically increases reaching, for all formulations, values that are about two orders of magnitude smaller than C p . This is essentially due to the high value of tissues volumes (V T = 24000 cm 3 ). Interestingly, for smaller particles, C p falls inside the theophylline therapeutic window that corresponds to 5 -15 g/cm 3 .
The second drug considered, temazepam (TEM), a benzodiazepine derivative characterised by sedative and hypnotic actions, was chosen for the possibility of loading it into a polymeric carrier (such as crosslinked polyvinyl-pyrrolidone particles) in a completely amorphous form [11]. This, of course, can be highly advantageous as amorphous TEM is characterised by a much higher solubility (≈ 8649 g/cm 3 ) with respect to the original macro-crystal solubility (≈ 164 g/cm 3 ). Accordingly, we would like to theoretically evaluate the PK advantage of using amorphous TEM in place of its native macro-crystal condition. In absence of precise values of TEM PK parameters, we adopted those of a similar drug [10] while all other model parameters were set according to our previous study [11] (see Table 2). Figure 3A makes clear that amorphous TEM (thick line) shows, in PK terms, a better performance with respect to crystalline TEM (thin line). Indeed, not only it yields to a higher maximum C P (0.16 g/cm 3 amorphous TEM, 0.14 g/cm 3 , crystalline TEM) but, also, maximum C P occurs after about 2 h instead of about 3.5 h (crystalline TEM). This simulation was performed assuming that amorphous TEM re-crystallisation did not occur inside the stomach although we have evidence of the contrary in water [11]. Actually, as re-crystallisation kinetics can strongly depend also on the physical/chemical properties of the dissolving fluids (stomach fluids differ a lot from simple water), lacking precise information about TEM re-crystallisation in the stomach fluids, induced us to assume that re-crystallisation did not occur. In so doing, we estimated the maximum possible effect of amorphous TEM over the native crystalline TEM. Figure 3B shows that, in the amorphous case, drug 90 concentration (C r ) rapidly reaches the maximum admissible concentration (400 g/cm 3 ) compatible with the administered dose (100 mg). This, obviously, takes place in the stomach (t < 0.8 h) while, once in the SIL and, eventually, in the LIL, C r decreases up to disappearing. In the crystalline case, on the contrary, C r increase is smaller and the maximum value is well below that reached in the amorphous case.

8649
At the same time, C r reduction in the SIL and in the LIL is slower thus ensuring higher C p from 4 hours onwards. In the case of a system characterised by 50 % amorphous TEM and 50 % crystalline TEM, both the C p and the C r trends are, more or less, in between those competing to the completely amorphous and the totally crystalline TEM as depicted by the grey lines in Figure 3A and 3B. Finally, it is interesting to notice that C r fastest variations occur when the delivery system leaves the stomach (0.8 h) and the LIL (3.2 h) ( Figure 3B) and this is determined by the GI tract permeability variation.
The last drug considered is nimesulide, a typical non-steroidal anti-inflammatory drug that has been largely employed in clinical practice [30]. Like TEM, it is scarcely soluble in aqueous media but its solubility is pH dependent as it is a weak acid (pKa = 6.46 [31]). Thus, water solubility at 37 °C is around 10 g/cm 3 for pH < pKa, while solubility increases to about 100 g/cm 3 with higher pH [32]. Obviously, nimesulide solubility can be improved by using it in its amorphous/nano-crystalline state [19]. In addition, the firstpass liver effect for nimesulide is negligible [30]. While some PK parameters [10,16,30] and all the delivery parameters [11,19] can be deduced from literature, we could not find reliable estimation of its solubility in the gastric environment and its permeability through the GL, SIL and LIL mucosa. Thus, our model was fitted (four fitting parameters) to experimental data referring to the oral administration (three humans) of a commercial prompt action formulation containing nimesulide (100 mg) [33]. Fitting procedure was performed assuming that nimesulide is embodied in macro-crystal form inside a polymeric carrier showing the physical-chemical characteristics reported in Table 3 (they are similar to those competing to polyvinyl-pyrrolidone particles considered in the TEM case). In addition, i GL R was set to 7.3 cm so that the stomach absorption surface equals the one normally found in humans (around 1 m 2 [16]). Figure 4A reports the very good agreement between model best fitting (solid thin line) and experimental data (open circles) being fitting parameters mc s C = 39 g/cm 3 , P GL = 6*10 -4 cm/s, P SIL = 8.3*10 -5 cm/s and P LIL = 7*10 -5 cm/s. These permeability values indicate that nimesulide is rapidly absorbed by the gastric mucosa while its absorption in the SIL and LIL is slower due to the reduction of permeability (one order of magnitude). These results seem reasonable. Indeed, the acid environment of the stomach is particularly favorable for the absorption of nimesulide that, at low pH (< 3; nimesulide is a weak acid with pKa = 6.46) is completely un-dissociated (un-ionized) and this condition favors its passive absorption by the lipid mucosal membrane. On the contrary, for higher pH, such as that found in the small intestine, nimesulide is partially or totally dissociated (ionized) and this condition hinders nimesulide crossing of the lipophilic mucosal membrane. Accordingly, chemical and physical properties should be responsible for the nimesulide different permeability through the GI tract rather than other reasons such as the effect of intestinal wall transporters. These chemico-phyiscal considerations have been supported by studies about the regional absorption of nimesulide [34,35], led on healthy volunteers (gamma scintigraphy). These studies proved that stomach and proximal small bowel account for 40 % of nimesulide absorption.  As in the TEM case, nimesulide absorption can be improved by considering its amorphous form. Figure  4A shows the predicted plasma concentration C p (thick line) in the case of amorphous nimesulide assuming, as done in the TEM case, that no re-crystallisation takes place (K rb = 0). Thus, this prediction should represent the maximum expected effect of amorphous nimesulide. It can be seen that the maximum concentration is considerably increased and, up to about 5 hours, drug concentration (thick line) is higher than that of the macro-crystal case (thin line). This behavior in the plasma is explained by Figure 4B where drug concentration in the release environment (C r ) is reported for the amorphous drug (thick line) and the macro-crystalline drug (thin line). Basically, in the amorphous case, nimesulide is rapidly released in the first 6 minutes and its concentration becomes about one order of magnitude higher than that corresponding to the macro-crystal case. This means that its onset of action is considerably increased, a desirable feature for a rapid pain relief. Also in this case, the time evolutions of drug concentration in the CIGS (C GICS ) and in the liver (C L ) are substantially equal to the C p one and the reasons are the same as the theophylline case's.
Again, drug concentration (C T ) in the scarcely perfused tissues and organs monotonically increases reaching, for all formulations, values that are about two order of magnitude smaller than C p . The abrupt increase of C r occurring after 1 h, in both the crystalline and amorphous case, is simply due to the reduced value of drug permeability in the LIL. Indeed, permeability reduction implies the decrease of the drug flux doi: 10.5599/admet.2.2.34 93 leaving the release environment and this, in turn, determines a drug concentration increase in the release environment.

Conclusions
Although the problem of a reliable prediction of drug concentration in plasma following oral administration remains a challenging task, a considerable development in this direction is represented by mathematical models able to simultaneously account for the in vivo drug release and the subsequent absorption, distribution, metabolism and elimination processes (ADME). In this frame, we built up a reasonable mathematical model resulting from merging a physiologically oriented PK model and a delivery model. Wherever the final aim is to get a reliable in vivo simulation, at present, we set up a theoretical tool that enables to rationally compare either different formulations of the same drug or the same formulation for different drugs. Indeed, our model allows the evaluation of the effect of different doses, different mean particles size and particle sizes distribution, and the drug solid states, i.e. amorphous, nano-crystalline and macro-crystalline drug. According to us, our model demonstrates that in vivo release kinetics can be different from the in vitro one due to the effect of living tissues. In other words, while the mutual effect between drug release and ADME processes is well known, this paper represents one of the first attempts of clearly evaluating and quantifying it. In conclusion, our model represents a rational tool for the designing of modern delivery systems and it should constitute an important step on the way to reliable simulation of plasma drug concentration following oral administration.