Mathematical modelling of drug permeation through a swollen membrane

Journal of Controlled Release 59 (1999) 343–359

Mathematical modelling of drug permeation through a swollen membrane a, 1 ,b M. Grassi *, I. Colombo a

Department of Chemical, Environmental and Raw Materials Engineering – DICAMP, University of Trieste, Piazzale Europa 1, I-34127 Trieste, Italy b Vectorpharma S.p. A., Via del Follatoio 12, I-34148 Trieste, Italy Received 19 August 1998; received in revised form 2 December 1998; accepted 3 December 1998

Abstract This work proposes two different mathematical models (linear and numerical) able to simulate the drug permeation through a swollen membrane sandwiched by two external layers (trilaminate system). Moreover, a solid drug dissolution phenomenon in the donor compartment may be accounted for. Indeed, this is a situation that may often occur in permeation experiments. An insufficient stirring of the donor and of the receiver volume may give rise to two sandwiching layers and the target of a constant drug concentration in the donor compartment may be accomplished by putting a solid drug amount in the saturated donor solution. The linear model shows the advantage of having an analytical expression which extremely simplifies the calculation of the drug diffusion coefficient D inside the membrane. Its main drawback lies in the fact that it works only for thin trilaminate systems. The numerical model is more general than the linear one, as it works for all kind of trilaminate thickness and it may account for a solid powder dissolution in the donor compartment. Of course, it does not have an analytical solution and, thus, the D determination is less easy to perform as the numerical model is more time consuming than the linear one. These two models are then compared with the classical approach developed by Flynn and Barrie in order to better define its validity limits.  1999 Elsevier Science B.V. All rights reserved. Keywords: Permeation; Mathematical modelling; Diffusion; Dissolution; Powder

1. Introduction In order to design a controlled release system based on a swellable polymeric matrix, it is of paramount importance to know the drug diffusion coefficient D in the matrix and its dependence on *Corresponding author. Tel.: 139-40-6763-435; fax: 139-40569-823. E-mail address: [email protected] (M. Grassi) 1 Present address: Eurand International S.p.A., Via Martin Luther King 13, I-20060 Pessano con Bornago, Milan, Italy.

temperature and on polymer and drug concentration [1]. Thus, the measure of D plays a very important role and the need to develop a proper model able to interpret the experimental data arises. Several methods are available in the literature for the experimental determination of D [2–4]. Among this plethora, we can mention the category of methods deriving from nuclear magnetic resonance (NMR) and dynamic light scattering (DLS) experiments, those based on holographic relaxation spectroscopy, those founded on the determination of the drug concentration profile such as the sectioning and inverse

0168-3659 / 99 / $ – see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S0168-3659( 98 )00198-9

344

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

sectioning method [5] and the methods based on drug concentration gradient under stationary and non-stationary gradients. In this paper we wish to focus our attention on the mathematical modelling of the drug permeation through a swollen membrane with particular attention to the possibility that two additional layers may sandwich the membrane itself and that a drug dissolution phenomenon may take place in the donor compartment. Three different approaches to the problem will be shown and they will be compared each other in order to stress their advantages and drawbacks.

2. Modelling One of the most important errors affecting the determination of the value of the drug diffusion coefficient obtained from permeation data through a swollen membrane, may be due to the presence of two stagnant layers arising, respectively, between the membrane and the bulk of the donor and receiver phase because of insufficient stirring. Neglecting the two layers means to determine the value of the diffusion coefficient, referred to the whole trilaminate (made up by the two stagnant layers and the membrane) instead of that referred to the single membrane. This implies an error depending on the sum of the thickness of the two stagnant layers. Sometimes this error may be not negligible, as the membrane thickness may be small. For instance, Giovannini [6] determined, in a classical side by side permeation apparatus, that the thickness of each stagnant layer may be around 60 mm, despite the presence of the magnetic stirrers in both the donor and receiver compartments. This thickness may be comparable with that of some common membranes employed in permeation experiments [7–9]. The mathematical models here considered, represent three different ways of taking into account the presence of the layers. One of the targets of this paper is to analyse these models and to evaluate how they behave in response to the variation of some fundamental parameters ruling the drug permeation. Of course, the main assumption over which they are built is that Fick’s law for the drug flux of matter

holds inside the trilaminate. This means that the effects of all possible chemical or electrical interactions between drug molecules and polymer chains or solvent molecules are accounted for by only considering a drug partition coefficient (membrane / release environment medium) different from the value of one. Of course, a more detailed analysis should account for drug–polymer interaction by coupling, for instance, the diffusion with a drug adsorption– desorption phenomenon on polymer chains [10,11]. Anyway, the usual way to proceed [12–16] is to implicitly incorporate the polymer–drug interactions in the drug partition coefficient. We believe that, usually, this assumption is more than enough to correctly model a permeation experiment. Furthermore, we suppose that the membrane is completely swollen, which means that the solvent concentration has reached its thermodynamic equilibrium value before starting the permeation experiments. Indeed, the membrane swelling may heavily influence the features of the drug permeation [17]. Fig. 1 schematically shows the physical set-up which all the following considerations will be referred to.

Fig. 1. Physical set-up. A membrane is sandwiched among two layers arising in the donor and receiver compartments due to an insufficient stirring. In the donor volume, a dissolution process may take place.

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

2.1. Linear model The main assumption adopted to formulate the linear model relies on the hypothesis that the concentration profile inside the two stagnant layers and the membrane has always a linear trend and that the drug diffusion coefficient is concentration independent. To infer a more general validity to this model, we suppose that the value of the drug concentration in the donor compartment is ruled both by the drug matter flux entering in the first stagnant layer and by the dissolution of a solid drug mass M (see Fig. 1). With such a hypothesis, the system of equations representing the linear model is the following one:

S

dCd C21 dM D1 S Vd ]] 5 2 ] 2 ]] Cd K1d 2 ] dt dt h1 K21

S

dCr D3 S C22 Vr ] 5 ]] ] 2 Cr K3r dt h 3 K23

S

D

D

(1)

(2)

D

D1 S C21 D2 S ]] Cd K1d 2 ] 5 ]]sC21 2 C22d h1 K21 h2 X 5 h1

S

D1 S D3 S C22 ]] C 2 C22d 5 ]] ] 2 Cr K3r h 1 s 21 h 3 K23 X 5 h1 1 h2 dM ] 5 2 KtsCs 2 CddVd dt

(3)

D (4) (5)

where Cd , Cr , Vd , and Vr are the drug concentration and the volume of the donor and the receiver compartment, respectively, h 1 , h 2 and h 3 represent the thickness of the first stagnant layer, of the membrane and the thickness of the second stagnant layer, respectively, D1 , D2 and D3 are the drug diffusion coefficients in the first layer, in the membrane and in the second layer, respectively, S is the available area for the permeation, Kt is the drug dissolution constant, C21 and C22 are the drug concentration in the membrane in X5h 1 and X5 h 1 1h 2 , respectively, t is time, Cs is the drug solubility in the donor and receiver fluid and K1d , K21 , K23 and K3r are the partition coefficients defined by the following equations:

345

C1sX 5 0d C1` K1d 5 ]]] 5 ] Cd Cd`

(6)

C2sX 5 h 1d C21 C2` K21 5 ]]] 5 ]]] 5 ] C1sX 5 h 1d C1sX 5 h 1d C1`

(7)

C2sX 5 h 1 1 h 2d C22 C2` K23 5 ]]]]] 5 ]]]]] 5 ] C3sX 5 h 1 1 h 2d C3sX 5 h 1 1 h 2d C3` C3sX 5 h 1 1 h 2 1 h 3d C3` K3r 5 ]]]]]] 5 ] Cr Cr`

(8) (9)

with C1 , C2 and C3 , respectively, the drug concentration in the first layer, in the membrane and in the second layer at the generic time t, while Cd` , Cr` ,C1` , C2` , C3` are the drug concentrations in the donor and receiver volume and the drug concentration in the 1 st layer, in the membrane and in the second layer, respectively, after an infinite time. In this way, we admit that the partition coefficients, defined by Eqs. (6)–(9), are concentration independent and, thus, that they are time independent. Eq. (1) represents the drug mass balance made up on the donor compartment: the first right hand side term takes in account the dissolution, while the second represents the matter flux leaving the donor trough the first stagnant layer. Eq. (2) represents the drug mass balance made up on the receiver compartment: the right hand side term is the entering drug flux coming from the second layer. Eq. (3) imposes that the matter flux leaving the first stagnant layer is equal to that entering the membrane (X5h 1 ), while Eq. (4) imposes the equality of the matter flux leaving the membrane and entering the second stagnant layer (X5h 2 1h 1 ). Finally, Eq. (5) rules the drug dissolution in the donor compartment. Two are the initial conditions to be satisfied by the solution of the above system of equations are: Crst 5 0d 5 0

(10)

Cdst 5 0d 5 Cd0

(11)

where Cd0 is the initial concentration in the donor compartment. The analytical form of the linear model may be achieved by firstly getting C21 and C22 as functions

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

346

of Cr and Cd by means of Eqs. (3) and (4), and then by solving the system of differential equations made up by Eqs. (1) and (2), provided that the right hand side of Eq. (5) is put in Eq. (1). The function M(t) may be obtained by integrating Eq. (5) knowing the function Cd (t). The final expressions for Cd (t), Cr (t) and M(t) are the following: Cd (t) 5 A 1 1 A 2 e sm 1 td 1 A 3 e sm 2 td m 1 td

Cr (t) 5 B1 1 B2 e s

m 2 td

1 B3 e s

m 1 td

M(t) 5 M0 1 E1se s

(12) (13) m 2 td

2 1d 1 E2se s

2 1d

(14)

where M0 is the M(t) starting value while A 1 , A 2 , A 3 , m 1 , m 2 , B1 , B2 , B3 , E1 and E2 are parameters depending on the geometrical and physical characteristics of the system made up by the donor, the two stagnant layers, the membrane and the receiver (see Appendix A for the explicit expression of these parameters). We have to stress the fact that Eqs. (12)–(14) implicitly admit that the dissolution phenomenon lasts until Cd and Cr reach their equilibrium value. This means that the amount of the beginning solid drug has to be sufficiently high. If it were not the case, we have to take into consideration the two forms of the model obtained by firstly setting Kt ±0 and then Kt 50. In order to generalise the results, the model prediction will be shown in terms of the following dimensionless variables: Cd Cr 1 tD2 1 1 C d 5 ]; C r 5 ]; t 5 ] Cd0 Cd0 h 22

(15) 2

Vd Vr Kt h 2 1 1 1 V d 5 ]; V r 5 ]; K t 5 ]] Sh 2 Sh 2 D2

(16)

As shown before, many are the parameters affecting C d1 and C r1 , but, in this work, we are mainly interested in investigating the effects of K t1 and the ratio h 1 /h 2 , considering, for sake of simplicity, h 1 5 h3. Fig. 2 shows the model prediction obtained assuming K 1 t 50, (D 1 /D 2 )5(D3 /D2 )51.5, K1d 5 K3r 51, K21 5K23 50.8, V d1 5V r1 51000, Cd0 5Cs , and four different values of H5h 1 /h 2 5h 3 /h 2 . The effect of H on C r1 and C d1 is considerably high: at a

1 Fig. 2. Effect of the layer thickness on the C 1 time d and C r profiles. As H increases, the permeation phenomenon is shifted towards longer times. The simulation is led assuming K 1 t 50, (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, V d1 5V r1 5 1000 and Cd0 5Cs .

fixed value of t 1 , its increase determines a strong decrease of the permeated mass. The thinner the membrane, the more important the effect of the stagnant layers on drug permeation because of the increase of H. The fact that C r1 and C d1 approach 0.5 is due to the fact that V d1 5V r1 , the volume of the trilaminate barrier being negligible. Fig. 3 reports the same situation depicted in Fig. 2 24 except for the fact that K 1 and, for a t 51.88?10 better understanding of the picture, only three values of H are reported. The effect of K t1 is clear: C r1 increases until reaching the final value of 1 because the drug amount contained in the donor, at the beginning of the permeation, is sufficient to saturate, at equilibrium, both the donor and the receiver. Initially, the matter flux leaving the donor through the trilaminate is higher than that coming in from the dissolution of the solid drug. So, an initial decrease

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

Fig. 3. Coupled effect of the dissolution and the layer thickness variation. After a starting sharp decrease, C 1 approaches the d solubility threshold, while C 1 increases monotonically towards r the same threshold. The simulation is led assuming K t1 51.88? 10 24 , (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, 1 V1 d 5V r 51000 and C d0 5C s .

of C 1 d takes place. The decrease of the concentration gradient across the trilaminate (due to a C 1 r increase in the receiver), reduces the matter flux leaving the donor until it becomes smaller than the dissolution one. In this way, an increase of C 1 d takes place until reaching the saturation value. The higher the K 1 t , the smaller the initial decrease and the faster the sub1 sequent increase of C 1 d . For very high K t values, 1 the minimum of C d curve disappears and a constant C1 d curve occurs. Also in this situation, the effect of 1 the ratio H on the C 1 r and C d trend is considerable: its increase determines a reduction of the matter flux through the trilaminate and, consequently, a lowering 1 of C 1 r at fixed t . Fig. 4 shows the drug concentration profile inside the trilaminate corresponding to three different values of t 1 , setting H50.5 with all other parameters

347

Fig. 4. Concentration profile inside the trilaminate system for three different t 1 values, assuming the same conditions set in Fig. 2.

equal to those of Fig. 2. The thinner line represents the profile concentration arising after a long time (t 1 55333) (thermodynamic equilibrium concentration), the thicker line indicates the profile concentration for t 1 50 and the middle line corresponds to the profile concentration developing for t 1 5533. It is interesting to notice that the profile concentration shows a discontinuity, for every t 1 , at the two membrane interfaces because of setting K21 5K23 5 0.8. For t 1 50, the model yields to a zero value for C r1 even if the profile concentration is completely developed. This is not reasonable since the drug needs some time to fill the membrane and, as a consequence, the gradient concentration does not develop instantaneously. This fact comes out directly from the hypothesis on which the model has been built up. For long permeation times, and thin trilaminate system, this is not a serious drawback because the filling time will be negligible but, for short permeation times, it may give origin to some prob-

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

348

lems. A possible way to empirically overcome this problem is to introduce a lag time t r taking in account the filling time. Eqs. (12) and (13) then become: Cd (t) 5 A 1 1 A 2 e sm 1st 2t rdd 1 A 3 e sm 2st2t rdd m 1st 2t rdd

Cr (t) 5 B1 1 B2 e s

m 2st 2t rdd

1 B3 e s

(17) (18)

Of course, Eqs. (17) and (18) may be employed only when we need to analyse an experimental data set related to a permeation time very small compared with that required to approach the equilibrium concentration both in the donor and in the receiver. The effect of t r is to cause a forward shifting of the C r1 curve vs. time. As an example, we would like to consider some experimental data determined by Giovannini [6] in a classical side by side permeation experiment. Briefly, she studied the permeation of theophylline through sodium alginate membranes of different thickness. The donor compartment contained a water saturated theophylline solution in presence of a solid, not dissolved, theophylline amount, while the receiver was initially filled with pure water. Fig. 5, showing the best fit of Eq. (18) (solid line) on a set of the above mentioned experimental data (open circles), underlines the good agreement between the model and the experiment. It is important to notice that, in this case, the H ratio (5h 1 /h 2 5h 3 /h 2 ) is equal to 0.21, this indicates a considerable weight of the two stagnant layers. On the contrary, Fig. 6, showing the best fit of Eq. (18) (solid line) on a different set (thicker membrane) of the above mentioned data (open circles), considers a situation where the effect of the stagnant layers is not so important, H50.09. The two different H values were obtained by increasing the membrane thickness, the stagnant layer thickness being constant, under the same hydrodynamic conditions for the two sets considered. These two examples demonstrate how the linear model is able to describe some different experimental conditions. It is now important to stress the advantages and the drawbacks connected to linear model use. As far as the advantages are concerned, it has to be said that the model has an analytical expression, it does not require the sink conditions in the receiver and a

Fig. 5. Best fit of Eq. (18) (solid line) on the experimental data (open circles) regarding the theophylline permeation through a thin swollen sodium alginate membrane.

constant value of Cd and, finally, it is able to take into account the effects of a possible dissolution phenomenon developing in the donor compartment during permeation. The drawbacks are related to the thickness of the trilaminate system made up by the two stagnant layers and the membrane. Indeed, the model gives a correct simulation of the permeation only for thin trilaminate systems, for which the drug concentration profile inside each stratum reaches a linear profile very quickly and, finally, the amount of the drug contained in the trilaminate is negligible (the thinner the trilaminate, the more accurate the model prediction).

2.2. Numerical model This approach to permeation is the more general one, as it does not require any particular assumption regarding the profile concentration inside the trilaminate, the diffusion coefficient and the value taken by 1 C1 r and C d during permeation itself. Unfortunately,

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

349

≠Cd dM ≠C1 Vd ]] 5 ] 1 D1 S ]]u X 50 ≠t dt ≠X

(22)

dM ] 5 2Vd KtsCs 2 Cdd dt

(23)

≠C1 ≠C2 D1 ]]u X5h 1 5 D2 ]]u X5h 1 ≠X ≠X

(24)

≠C2 ≠C3 D2 ]]u X5h 1 1h 2 5 D3 ]]u X5h 1 1h 2 ≠X ≠X

(25)

≠Cr ≠C3 Vr ] 5 D3 S ]]u X5h 1 1h 2 1h 3 ≠t ≠X

(26)

C1sX 5 0d 5 K1d Cd

(27)

C2sX 5 h 1d 5 K21 C1sX 5 h 1d

(28)

C2sX 5 h 1 1 h 2d 5 K23 C3sX 5 h 1 1 h 2d

(29)

C3sX 5 h 1 1 h 2 1 h 3d 5 K3r Cr

(30)

and the following initial conditions: Cd 5 C1 5 Cd0 ; Cr 5 C2 5 C3 5 0; M 5 M0 Fig. 6. Best fit of Eq. (18) (solid line) on the experimental data (open circles) regarding the theophylline permeation through a thick swollen sodium alginate membrane.

this approach does not lead to an analytical solution making the model difficult to use. For these reasons, it has to be employed as a reference for other models and when the other models do not work. The numerical model is nothing more than the solution, in one dimension, of Fick’s second law inside the trilaminate undergoing the pertinent initial and boundary conditions. Fick’s second law for the first layer, the membrane and the second layer reads, respectively:

S D ≠C ≠C ≠ ]] 5 ]SD ]]D ≠t ≠X ≠X ≠C ≠C ≠ ]] 5 ]SD ]]D ≠t ≠X ≠X ≠C ≠C ≠ ]]1 5 ] D1 ]]1 ≠t ≠X ≠X 2

2

(20)

3

(21)

2

3

3

(19)

where X is the abscissa. The above equations must be solved with the following boundary conditions:

(31)

Eq. (22) represents the drug mass balance made up on the donor compartment: the first right hand side term takes into account the dissolution, while the second represents the matter flux leaving the donor trough the first stagnant layer. Eq. (23) takes into account the reduction of the solid drug M as the dissolution goes on. When M is zeroed, the first term of the right hand side of Eq. (22) vanishes. Eq. (24) imposes that the matter flux leaving the first stagnant layer is equal to that entering the membrane (X5h 1 ), while Eq. (25) imposes the equality of the matter flux leaving the membrane and entering the second stagnant layer (X5h 2 1h 1 ). Eq. (26) represents the drug mass balance made up on the receiver compartment: the right hand side term is the entering drug flux coming from the second layer. Eqs. (27)–(30) indicate the partitioning relation holding at the interface in X50, h 1 , (h 1 1h 2 ), (h 1 1h 2 1h 3 ). Eq. (31) sets to zero the drug concentration in the membrane, in the second layer and in the receiver while it sets to Cd0 the drug concentration in the first stagnant layer and in the donor at the beginning of the permeation. As said before, the solution of the above set of

350

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

equations may be achieved only by means of a numerical method. We chose the control volume method [18] which is an implicit finite differences method suitable to solve such kind of problems [19–21]. It requires one to subdivide the three strata of the trilaminate in N1 , N2 and N3 elementary volumes, respectively, and to integrate, in the space and over time, Fick’s second equation. This model resembles that of Kurnik and Potts [22], even if it deals with a more complex geometry and it describes a little bit different physical situation. Indeed, while Kurnik takes into account a dissolution process developing inside the membrane, this model considers a dissolution process taking place only in the donor compartment. Fig. 7 shows the simulation results achieved by means of the numerical model assuming the same values of the parameters adopted in Fig. 2, setting the dimensionless time increment Dt 1 51.3325 and N1 5N2 5N3 530. The numerical model and the

linear one give the same prediction for both C 1 r and 1 1 1 C d since the parameter V r (5V d ) is sufficiently high, indeed, for a constant value of Vr and S, V r1 increases as h 2 decreases. So, high values of V 1 r correspond to a sufficiently small trilaminate thickness and a linear concentration profile in each part of the trilaminate may be assumed. Fig. 8 stresses this last consideration showing that also the numerical model yields to a linear concentration profile in each part of the trilaminate and, of course, this trend is equal to that of Fig. 4. Nevertheless, for very short time, the numerical model furnishes a clearly non linear concentration profile inside the trilaminate indicating the lack of physical description of the linear model at the first beginning of the permeation. The numerical model and the linear one yield to 24 the same prediction in the case of K 1 . t 51.88?10 It is easy to demonstrate that the lower V r1 (5V 1 d ), the bigger the differences in the prediction obtained from the linear and numerical model. These

Fig. 7. Numerical model simulation of the permeation through a trilaminate system adopting the same conditions set in Fig. 2. As 1 V1 are sufficiently big, the linear and the numerical d and V r model yield the same result.

Fig. 8. Concentration profile inside the trilaminate system for four different t 1 values, assuming the same conditions set in Fig. 2. It is clearly evident that, for the smallest time (t 1 50.266), the concentration profile is not linear.

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

differences are more evident at the beginning of the permeation since, for long time, a linear profile concentration will be completely developed in the three part of the trilaminate. Fig. 9 shows the two model beginning predictions for two different V 1 r 1 (5V 1 d ) values and assuming K t 50, (D 1 /D 2 )5(D 3 / D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, Cd0 5Cs and H50.5. It is important to notice that for V r1 5100 the two curves are parallel while, for V 1 r 520, they have the tendency to diverge. This is due to the fact that the drug amount contained in the trilaminate is no more negligible and the linear model, neglecting this amount, yields to an overestimation also of the C1 r slope curve. In this situation, indeed, the equilibrium value predicted by the two models is not the same. Considering the whole permeation time, assuming V1 r 5100 and the same parameters adopted in Fig. 9,

Fig. 9. Comparison between the numerical and linear model 1 prediction for two different values of V 1 r . For V r 520, the two model simulations have the tendency to diverge as a consequence of the higher importance of the mass accumulation in the trilaminate. The simulation is led assuming K 1 t 50, (D 1 /D 2 )5 (D3 /D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, H50.5 and Cd0 5Cs .

351

the linear and numerical models yield to an almost similar prediction, being the filling time of the trilaminate and the drug amount contained in it small. This can be seen in Fig. 10 with the percentage relative difference (RD) defined by the following equation: 1

C r,lin 2 Cr,num % RD 5 100]]]] Cr,num

(32)

where the subscripts ‘‘lin’’ and ‘‘num’’ refer to the C1 value obtained by the linear and numerical r model, respectively. It is interesting to notice that RD is always lower than 9% and, for t 1 greater than 60, it is lower than 2%. Up until now we supposed to know the value of K1 although its knowledge may require many t detailed information about the dissolution phenomenon. Indeed, the dissolution of a solid from a plane and uniform surface may be easily modelled by means of the following equation [23]:

Fig. 10. Percentage relative difference (RD) between the linear and the numerical model simulation as a function of t 1 . It is assumed that K 1 (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, t 50, 1 K21 5K23 50.8, H50.5 and Cd0 5Cs and V 1 r 5V d 5100.

352

dC D0 Sp V ] 5 ]]sCs 2 Cd dt h

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

(33)

where C, Cs and V are, respectively, the solute concentration, the solute solubility and the volume of the dissolution medium, Sp is the area of the solid– liquid interface, D0 is the solute diffusion coefficient in the dissolution medium and h is the thickness of the boundary layer arising between the solid surface and the dissolution medium. h strongly depends on the stirring conditions of the dissolution medium as stated by Levich [24] and Banakar [25]. Remembering that Kt is equal to: D0 Sp D0 K1 5 ]]; Kd 5 ] hV h

(34)

The use of the numerical model allows one to also consider the above mentioned case of a powder dissolution characterised by a time dependent Kt . In order to match this purpose, we have to insert the Kt time dependency, given by Eq. (38), into Eq. (22), getting: ≠Cd ≠C1 Vd ]] 5 4p Np R 2 KdsCs 2 Cdd 1 D1 S ]]u X50 ≠t ≠X

This equation has to be coupled with the R time dependency coming out from the following particle mass balance:

S

(35)

The right hand side of Eq. (35) coincides with the dissolution contribute employed in Eq. (22). Since Sp , for a plane and uniform surface, does not change as dissolution develops, we may be sure that Kt is time independent. Unfortunately, this is not the case for a dissolving powder. Indeed, in such hypotheses, the dissolution surface decreases as the time goes on. The starting value of Sp , Sp0 , is equal to the powder surface area and, for a monodisperse powder made up by Np , all equal, spherical particles, is given by: Sp0 5 Np 4p R 02

(36)

where R 0 is the particle radius. At time t, the particle radius will be decreased to R and, as a consequence, Sp will be equal to Sp 5 Np 4p R 2

(37)

Then, the Kt time dependency will be given by: Kd 2 K t 5 ] 4p R V

(38)

Eq. (38) holds in the hypothesis that the whole powder surface is available for dissolution. This is reasonably accomplished when the dissolution medium is highly stirred so that the particles cannot adhere to the vessel walls or each other.

D

3 dMp d 4pr R dR ]] 5 ] ]] 5 4pr R 2 ] 5 dt dt 3 dt

Kd 4p R 2sCs 2 Cdd

where Kd is the drug dissolution rate. Eq. (33) may be recast in the following form: dC V ] 5 KtsCs 2 CdV dt

(39)

(40)

where r is the particle density and Mp is the mass particle at time t. From Eq. (40) we get [24]: dR K ] 5 ]d sCs 2 Cdd dt r

(41)

modelling the R time dependency. Despite the theoretical correctness of Eq. (41), because of numerical problems, we prefer to evaluate the R variation by making recourse to a mass balance made up on the trilaminate, the donor and receiver compartments. This mass balance reads: h1

E

M 5 M0 1VdsCd0 2 Cdd 2Vr Cr 2 C1 SdX 0 h 1 1h 2

2

h 1 1h 2 1h 3

E C SdX 2 E 2

h1

C3 SdX

(42)

h 1 1h 2

with C1 , C2 , C3 , h 1 , h 2 and h 3 , respectively, the drug concentration and the thickness of the first layer, the membrane and the second layer of the trilaminate system. Bearing in mind that: 4 M 5 Np Mp 5 Np ]pr R 3 3 we have: ]] 3M 3 R 5 ]] Np 4pr

œ

(43)

(44)

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

353

The use of Eq. (44) allows to assume a wider time step in the numerical solution [18] of our set of equations (Eqs. (23)–(31), (39), (44)) avoiding stability and accuracy problems. An easier expression of Eq. (44) may be obtained remembering that, for a monodisperse powder made up by N spherical particles, the following relations must hold: 4 M0 5 Np ]p R 30 r ; M0 A 5 Np 4p R 20 3

(45)

where A is the powder surface per unit mass. The solution of the above system leads to: 3M0 3 R 0 5 ]; Np 5 ]] Ar 4p R 03 r

(46)

By using Eq. (46), Eq. (44) may be recast in the following form: ] M 3 R 5 R0 ] (47) M0

œ

which is the well known Hixon Crowell equation [26]. The effects of the Kt reduction are more evident when the dissolution phenomenon implies a considerable decreasing of the particle radius. Indeed, in this case, the dissolution surface will be strongly reduced and, as a consequence, the dissolving drug mass going into the donor compartment will be decreased. Fig. 11 shows the simulation results obtained supposing to neglect the K 1 reduction t (solid lines) and supposing to consider the K t1 1 24 reduction (dashed lines). We set K t 51.88?10 , (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, K21 5K23 5 0.8, H51, Cd0 5Cs , V d1 51000, V r1 58000, M 01 5 0.8, r 1 5120 and A1 5400.2, where: M0 Ah 2 r 1 1 1 M 0 5 ]]; r 5 ]; A 5 ]; Cd0Vd Cd0 Cd0 A g 5 AM0

(48)

1 For each case, both C 1 d (the upper line) and C r (the lower line) are depicted. It is clearly evident how taking into account the K 1 t reduction reflects in 1 a slower C 1 r and C d increase after the initial sharp reduction. This considerable difference is essentially due to the fact that, having set a great value for V 1 r , the dissolution phenomenon requires a great reduc-

1 1 Fig. 11. Effect of the K 1 t time reduction on C d and C r supposing 1 24 to set K t 51.88?10 , (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, h 1 5h 2 5h 3 , Cd0 5Cs , V d1 51000, V r1 58000, 1 1 M1 0 50.8, r 5120, A 5400.2.

tion of the dissolution surface reflecting in an equal reduction of K 1 t . We would have not had remarkable 1 1 differences in the C 1 d and C r curve if V r would have been equal to 1000: in this case the reduction of the dissolution surface would have been negligible. 1 Fig. 12 reports the K t reduction occurring during the permeation phenomenon simulated in Fig. 11 1 jointly with the K t reduction related to a lower A g value (1 / 4 reduction), keeping constant Kd . It results that, all other parameters being equal, the bigger A, the bigger the initial K 1 t value and the stronger the subsequent reduction. This tells us that the effect of the K t1 reduction becomes smaller and smaller as the initial powder surface per unit mass is lower. The numerical model has the advantage of considering the non-sink receptor conditions and the C 1 d variations during permeation jointly with high membrane thickness, drug diffusion coefficient dependence on drug concentration and dissolution from a powdered solid drug. Its main drawback lies in the

354

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

X C 5 KC0 ] hm

O

S D

`

KC0 2 np X 2((n 2 1 ] ]]cos (np )sin ]] e p / h m ) Dt ) p n51 n hm (49) where X is the one dimensional co-ordinate, h m is membrane thickness, K is the partition coefficient and t is time. In order to get Cr at time t, Eq. (49) has to be differentiated with respect to X to obtain the instantaneous drug flux in X5h m , then the drug flux has to be integrated from t50 to t * 5t and the result has to be divided for Vr . The resulting equation is: C0 KSDt C0 KSh m Cr 5 ]]] 2 ]]] hVr 6Vr

O

2h m C0 KS ` (21)n 2((n p / h m ) 2 Dt ) ]]e 1 ]]] p 2Vr n51 n 2

1 t

Fig. 12. Comparison between the K reduction occurring for two different values of the global surface area A g . The solid line represents the K 1 t reduction occurring in the simulation of Fig. 11. The dotted line corresponds to the K 1 t reduction occurring in the hypothesis of a global surface area equal to one-fourth of that fixed in Fig. 11.

fact that it has not an analytical expression, being not easy to use and being high time consuming in experimental data fitting. For these reasons, it is convenient to employ it when the linear and steady state model fail.

2.3. Steady-state model The well detailed analysis of the permeation phenomenon led by Flynn et al. [14] and Barrie et al. [15] allows one to formulate a steady-state model for a trilaminate system. Assuming the sink conditions in the receptor compartment, a constant drug concentration C0 in the donor compartment and a constant value of the drug diffusion coefficient D, the solution, in one dimension, of Fick’s second law, relative to a membrane, reads:

(50)

where S is the diffusing surface area. For t approaching infinity, Eq. (50) coincides with a straight line given by:

S

2 C0 KSD hm Cr 5 ]] t 2 ] h mVr 6D

D

(51)

Defining the membrane permeability P and the lag time t L as follows: h 2m KD P 5 ]; t L 5 ] hm 6D

(52)

Eq. (51) becomes: C0 SP Cr 5 ]]st 2 t Ld Vr where (51). For Flynn holds,

(53)

t L represents the time axis intercept of Eq. a trilaminate system, Barrie et al. [15] and et al. [14] demonstrated that Eq. (53) still but with different expressions for P and t L :

D1 D2 D3 K1 K2 K3 P 5 ]]]]]]]]]]]]] h 1 D2 D3 K2 K3 1 h 2 D1 D3 K1 K3 1 h 3 D1 D2 K1 K2 (54)

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359 tL 5

3 F

S

h 21 h1 h2 h3 ] ]] 1 ]] 1 ]] D1 6D1 K1 2D2 K2 2D3 K3 h 23

S

D S

h 22 h1 h2 h3 1 ] ]] 1 ]] 1 ]] D2 2D1 K1 6D2 K2 2D3 K3

h1 h2 h3 ] ]] 1 ]] 1 ]] D3 2D1 K1 2D2 K2 6D3 K3

h1 h2 h3 ]1]1] D1 K1 D2 K2 D3 K3

D

K2 h 1 h 2 h 3 1 ]] D1 D3 K1 K3

G

D

1

4

/

(55)

where K1 5 K1d ; K2 5 K21 K1d ; K3 5 K21 K3r /K23

(56)

Eqs. (54)–(56) allow to simulate the drug permeation phenomenon for a trilaminate at the steady state when the concentration gradient across the trilaminate is constant with time. Fig. 13 shows the comparison between the simulation results coming from the steady state 1 model and the numerical one assuming V 1 r (5V d )5 20, (D1 /D2 )5(D3 /D2 )51.5, K1d 5K3r 51, K21 5 K23 50.8, Cd0 5Cs , H50.5 and two different K 1 t

355

1 values. For K 1 decreases t 50 (dotted line), C d during permeation invalidating, in this manner, the basic assumption upon which the steady state model is made up. As a consequence, the steady state model and the numerical one disagree. Assuming a higher K1 value (K t1 510 2 solid thick line) serves to t increase the drug dissolution flux and, thus, it serves to strongly reduce the C 1 d decrease. In this new situation, indeed, the numerical model prediction approaches the steady state one, although the curve slope is still different. This is due to the fact that the sink conditions in the receiver are not accomplished because of the small value of V 1 r . This last consideration stress the fact that the steady state model suitability is very difficult to be matched when working with little receiver volume. Among the shown models, the steady-state one represents the simplest way of simulating, at stationary, the drug permeation phenomenon through a membrane (regardless of its thickness) sandwiched between two layers provoked by insufficient donor and receptor stirring. Anyway, as said before, it requires some assumptions that are not always easy to guarantee.

3. Conclusions

Fig. 13. Comparison between the steady state model prediction and the numerical one in the hypothesis of setting (D1 /D2 )5(D3 / D2 )51.5, K1d 5K3r 51, K21 5K23 50.8, h 1 /h 2 5h 2 /h 3 50.5, V r1 5 1 V1 d 520 and C d0 5C s . Two different values of K t are considered.

The above analysis demonstrates the necessity of verifying the presence and the entity of the stagnant layers sandwiching the membrane before measuring the drug diffusion coefficient from experimental permeation data. At the same time, it stresses the importance of an exact knowledge of Kt in order to properly evaluate the effect of the dissolution process on the drug concentration in the donor and in the receptor compartment. Indeed, with a sufficiently high Kt value, it is possible to have a constant value of Cd even in the presence of a relative small donor volume. In this situation a constant gradient concentration across the membrane arises, provided that the receptor volume is sufficiently big. This is the typical situation in which the steady-state model becomes the better tool in order to measure D because, being all the constrains it requires matched, it results to be faster and easier to adopt than the other two models. When Kt is not high enough and the matter

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

356

increase in the receptor compartment has to be taken into account, for thin membranes, the linear model becomes the most appropriate tool to measure D. It is also possible to adopt the linear model for thicker membranes by introducing a lag time t r in the time model dependence. Of course, the thickness must not to be so great as to make the mass of drug contained in the membrane with respect to the receiving one at the beginning of the permeation experiment negligible. The numerical model represents the only available tool when the other two models fail because of the lack of some fundamental physical requirements. This model has to be also used to check the reliability of the other model prediction.

C0 C1 C2 C3 C21

4. List of symbols

D1 D2 D3

a A A1 Ag A1 A2 A3 b B1 B2 B3 C Cd C d1 Cd0 Cd` Cr Cr` C r1 C1 r,lin 1 C r,num

Cs

Parameters (Eq. (59)) Powder surface per unit mass Dimensionless powder surface per unit mass Global powder surface Parameters (Eq. (69)) Parameters (Eq. (70)) Parameters (Eq. (71)) Parameters (Eq. (59)) Parameters (Eq. (72)) Parameters (Eq. (73)) Parameters (Eq. (74)) Drug profile concentration inside the membrane (Eq. (49)) Donor drug concentration Dimensionless donor drug concentration Starting donor drug concentration Donor drug concentration after an infinite time Receiver drug concentration Receiver drug concentration after an infinite time Dimensionless receiver drug concentration C1 r value calculated by means of the linear model C1 value calculated by means of the r numerical model Drug solubility in the donor and receiver fluid

C22 C1` C2` C3` D D0

E1 E2 g G h

H hm h1 h2 h3 K Kd Kt K t1 K1 K2 K3 K1d K21 K23 K3r M Mp M0 M 01

Constant drug concentration (Eq. (49)) Drug concentration in the first layer Drug concentration in the membrane Drug concentration in the second layer Drug concentration in the membrane in X5h 1 Drug concentration in the membrane in X5h 1 1h 2 Drug concentration in the first layer after an infinite time Drug concentration in the membrane after an infinite time drug concentration in the second layer after an infinite time drug diffusion coefficient (Eq. (50)) solute diffusion coefficient in the dissolution medium Drug diffusion coefficient in the first layer Drug diffusion coefficient in the membrane Drug diffusion coefficient in the second layer Parameter (Eq. (75)) Parameter (Eq. (76)) Parameter (Eq. (59)) Parameter (Eq. (63)) Thickness of the boundary layer arising between solid surface and dissolution medium h 1 /h 2 5h 3 /h 2 ratio Membrane thickness First layer thickness Membrane thickness Second layer thickness Partition coefficient Drug transfer rate Drug dissolution constant Dimensionless drug dissolution constant Partition coefficient (Eq. (56)) Partition coefficient (Eq. (56)) Partition coefficient (Eq. (56)) Partition coefficient (Eq. (6)) Partition coefficient (Eq. (7)) Partition coefficient (Eq. (8)) Partition coefficient (Eq. (9)) Solid drug mass Particle mass Starting solid drug mass Dimensionless starting solid drug mass

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

m1 m2 Np N1 N2 N3 P R R0 RD S Sp Sp0 t T 1 t tL tr V Vd V d1 Vr V1 r X x Y Z

Parameter (Eq. (67)) Parameter (Eq. (68)) Number of particles Number of elementary volumes in which the first layer is subdivided Number of elementary volumes in which the membrane is subdivided Number of elementary volumes in which the second layer is subdivided Membrane permeability Particle radius Starting particle radius Relative difference (%) Available area for permeation Solid–liquid interface area Starting solid–liquid interface area Time Parameter (Eq. (63)) Dimensionless time Lag time (Eqs. (52) and (55)) Lag time (Eqs. (17) and (18)) Dissolution medium volume Donor volume Dimensionless donor volume Receiver volume Dimensionless receiver volume Abscissa Parameter (Eq. (63)) Parameter (Eq. (63)) Parameter (Eq. (63))

4.1. Greek letters a b g d r r1 Dt 1

Parameter (Eq. (60)) Parameter (Eq. (60)) Parameter (Eq. (60)) Parameter (Eq. (60)) Particle density Dimensionless particle density Dimensionless time interval

357

and C22 on Cd and Cr . So, solving the linear system made up by Eqs. (3) and (4), one gets: C22 5 bCr 1 aCd

(57)

C21 5 zCr 1 gCd

(58)

where:

g b ag a 5 ]]; b 5 ]]; g 5 ]]; a 2d a 2d a 2d db z 5 ]] a 2d

(59)

and: D3 h 2 D3 h 2 K3r a 5 1 1 ]]]; b 5 2 ]]]; D2 h 2 K23 D2 h 3 K1d 1 g 5 ]]]]; d 5 ]]]] D2 h 2 D h 1 1 2 ]] 1 ] ]]] 11 D1 h 2 K21 D2 h 1 K21

(60)

Putting Eqs. (57) and (58) into Eqs. (1) and (2), one obtains the following system of differential equations: dCd ]] 5 Kt Cs 1 CdsTg 2 Gd 1 TzCr dt

(61)

dC ]r 5 Crsxb 2 Yd 1 xaCd dt

(62)

where: D1 SK1d D1 S G 5 Kt Cs 1 ]]; T 5 ]]]; Vd h 1 Vd h 1 K21 D3 S D3 S x 5 ]]]; Y 5 ]] h 3Vr K23 h 3Vr

(63)

The solution of the above system may be achieved by the usual technique suitable for system of linear first-order differential equations [27]. The solutions are: Cd (t) 5 A 1 1 A 2 e sm 1 td 1 A 3 e sm 2 td

(64)

Cr (t) 5 B1 1 B2 e sm 1 td 1 B3 e sm 2 td

(65)

M(t) 5 M0 1 E1se sm 1 td 2 1d 1 E2se sm 2 td 2 1d

(66)

Appendix A Some details about the solution of the system made up by Eqs. (1)–(5) are reported below. As a first step, we have to determine, by means of Eqs. (3) and (4), the functionality dependence of C21

where:

358

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359

m 1 5 0.5 f 2sTg 2 G 1 xb 2 Yd ]]]]]]]]]]]]] 1œsTg 2 Gxb 2 Yd 2 2 4 hsxb 2 YdsTg 2 gd 2 Tzxa j g

[5]

(67) m 2 5 0.5 f 2sTg 2 G 1 xb 2 Yd ]]]]]]]]]]]]] 1œsTg 2 Gxb 2 Yd 2 2 4 hsxb 2 YdsTg 2 gd 2 Tzxa j g (68) sxb 2 YdKt Cs A 1 5 2 ]]]] m1m2

[6]

[7]

(69)

sxb 2 YdKt Cs Kt Cs 2 ]]]] 2 f m 2 2sTg 2 Gd g Cs m1 A 2 5 ]]]]]]]]]]]] m1 2 m2

[8]

[9]

(70) sxb 2 YdKt Cs Kt Cs 2 ]]]] 2 f m 1 2sTg 2 Gd g Cs m2 A 3 5 ]]]]]]]]]]]] m2 2 m1

[10]

(71)

[11]

Kt Cs 1sTg 2 Gd A 2 B1 5 2 ]]]]]] Tz

(72)

A2 B2 5 ] f m 1 2sTg 2 Gd g Tz

(73)

[13]

A3 B3 5 ] f m 1 2sTg 2 Gd g Tz

(74)

[14]

A2 E1 5 ] m1

(75)

A3 E2 5 ] m2

(76)

[12]

[15]

[16]

[17]

References [18] [1] M. Tanya am Ende, N.A. Peppas, Investigation of drug / polymer interactions in the hydratated state, Proceedings Conference on Advances in Controlled Delivery, Baltimore, MD, 19–20 August 1996, pp. 12–13. [2] B.A. Westrin, A. Axelsson, G. Zacchi, Diffusion measurement in gels, J. Control. Release. 30 (1994) 189. ¨ [3] L. Johansson, J.E. Lofroth, Diffusion and interactions in gels and solution. I. Method, J. Colloid Int. Sci. 142 (1991) 116. [4] S.K. Inoue, R.B. Guenther, S.W. Hoag, Algorithm to determine diffusion and mass transfer coefficients, Proceedings

[19] [20]

[21] [22]

Conference on Advances in Controlled Delivery, Baltimore, MD, 19–20 August 1996, pp. 145–146. I. Colombo, M. Grassi, R. Lapasin, S. Pricl, Determination of drug diffusion coefficient in swollen hydrogel polymeric matrices by means of the inverse sectioning method, J. Control. Release 47 (1997) 305. F. Giovannini, Sistemi farmaceutici a rilascio controllato basati su gel di alginato, Graduate Thesis, Universita` di Trieste, Facolta` di Ingegneria, Dipartimento di Ingegneria Chimica (DICAMP), 1996. Y. Maiatami, H. Sato, T. Nagai, New method of calculating the diffusion coefficient of drugs from permeability through a membrane using diclofenac and its sodium salt in ethanol aqueous solution, Proc. 21st Int. Symp. on Control. Release of Bioactive Materials, 1994, p. 1234. G. Imanidis, S. Helbing-Strausak, R. Imboden, H. Leuenberger, Vehicle-dependent in situ modification of membranecontrolled drug release, J. Control. Release 51 (1998) 23–34. J.R. Nixon, K.T. Wong, Evaluation of drug permeation through polymeric membrane as a model for release (II) ethylcellulose-walled microcapsules, Int. J. Pharm. 58 (1990) 31–40. M. Singh, J.A. Lumpkin, J. Rosenblatt, Mathematical modeling of drug relese from hydrogel matrices via a diffusion coupled with desorption mechanism, J. Control. Release 32 (1994) 17. A. Adrover, M. Giona, M. Grassi, Controlled release of theophylline from water-swollen scleroglucan matrices, J. Membr. Sci. 113 (1996) 7. J. Crank, The Mathematics of Diffusion, Clarendon Press, Oxford, 1975. N.A. Peppas, in: R. Langer and D. Wise (Eds.), Medical Applications of Controlled Release: Applications and Evaluations, Vol. II, CRC Press, Boca Raton, FL, 1984. G.L. Flynn, S.H. Yalkowsky, T.J. Roseman, Mass transport phenomena and models: theoretical concepts, J. Pharm. Sci. 63(4) (1974) 479. J.A. Barrie, J.D. Levine, A.S. Michaels, P. Wong, Diffusion and solution of gases in composite rubber membranes, Trans. Faraday Soc. 59 (1962) 869. A. Axelsson, B. Westrin, Application of the diffusion cell for the measurement of diffusion in gels, Chem. Eng. Sci. 46(3) (1991) 913. M. Grassi, S.H. Yuk, S.H. Cho, Modelling of solute transport across a temperature-sensitive polymer membrane, J. Membr. Sci., in press. S.V., Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing, New York, 1980. G. Camera-Roda, G.C. Sarti, Mass transport with relaxation in polymers, AIChE J. 36(6) (1990) 851. M. Grassi, R. Lapasin, S. Pricl, Drug release from polymeric matrices with different degree of swelling, Proceedings of The First European Congress on Chemical Engineering, Firenze, 4–7 May 1997, pp. 2093–2095. M. Grassi, R. Lapasin, S. Pricl, Modelling of drug release from a swellable matrix, Chem. Eng. Commun., in press. R.T. Kurnik, R.O. Potts, Modeling of diffusion and crystal

M. Grassi, I. Colombo / Journal of Controlled Release 59 (1999) 343 – 359 dissolution in controlled release systems, J. Control. Release 45 (1997) 257–264. [23] H. Nogami, T. Nagai, E. Fukoka, T. Yotsuyanagi, Dissolution kinetics of barbital polymorphs, Chem. Pharm. Bull. 17(1) (1969) 23. [24] V.G. Levich, Physicochemical Hydrodynamics, PrenticeHall, Englewood Cliffs, NJ, 1962, pp. 39–73.

359

[25] U.V. Banakar, Pharmaceutical Dissolution Testing, Marcell Dekker, New York, Basel, Hong Kong, 1992, pp. 27–36. [26] U.V. Banakar, Pharmaceutical Dissolution Testing, Marcell Dekker, Inc., New York, Basel, Hong Kong, 1992, p. 10. [27] B. Demidovic, Esercizi e Problemi di Analisi Matematica, Editori Riuniti, Rome, 1975.