L J Fermín

J Lévy Véhel

Variability and singularity arising from poor compliance in a pharmacokinetic model II: the multi-oral case

 and help understanding the consequences of poor compliance. They may have practical outcomes in terms of drug dosing and scheduling.

Introduction

In the seminal work [START_REF] Li | A Pharmacokinetic Formalism Explicitly Integrating the Patient Drug Compliance[END_REF], the authors attacked the problem of mathematically modeling poor compliance using a probabilistic frame. They considered general distributions for the random instants of drug intake and studied the mean and variance of the concentration conditioned on the time elapsed since the last intake. Various related models were studied in detail in [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF]: the probability distribution of drug concentration in the context of multiple-IV dosing and poor compliance was investigated when both the times of intakes and doses are supposed to be random. In addition, both continuous-time and discrete-time versions were considered.

The present work parallels the studies of [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF] in the more realistic framework of multiple oral dosing. As in [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF], we investigate below the probability distribution of drug concentration in the context of poor compliance: the instants of drug intake are randoms, with possibly random doses. We use the simplest possible law to model the random times of drug intake, i.e. a homogeneous Poisson distribution. In other words, the times of drug intake are supposed to follows a Poisson process. This assumption allows one to perform explicit computations using the well-developed machinery on Poisson processes, and to obtain precise results describing various aspects of the distribution of the concentration that are important for assessing the efficacy of the regimen. We focus on two aspects of practical relevance: the variability of the concentration and the regularity of its probability distribution.

Our results quantify the variability of the concentration around its mean showing the exact role played by each parameter of the process. We measure how much poor compliance increases this variability as compared to the full compliance case. We show that the probability distribution limit of drug concentration may display irregular behaviours: in our case this occurs when the ratio between the mean number of intakes per unit time and the minimum between the elimination rate and the absorption rate is smaller than one half. When the same ratio is smaller than one, it becomes singular (i.e. non differentiable) at the origin. In practical terms, this amounts to quantifying, in a precise way, the situations where the moments of intakes are too scarce (with respect to the elimination and absorption rates), resulting in a high probability of having too small a concentration of drugs.

We also study a discrete-time version of the model. This setting reveals unexpected links with multifractal measures which have been studied in the mathematical literature for over seventy years under the name of infinite Bernoulli convolutions. Again, depending on some parameters, the discretized concentration may exhibit an extremely irregular behaviour. This means that the probability of observing a concentration C depends in a very non-smooth way on the precise value of C. This is obviously an undesirable feature which may have strongly negative consequences.

As mentioned above, the present article is the second in a series of three, where the first work [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF] dealt with the multi IV case. Since some computations here are very similar to this case, they are omitted and the interested reader is referred to [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF] for full details. The third paper in the series, [START_REF] Fermin | Modeling patient poor compliance in the multi-IV administration case with Piecewise Deterministic Markov Models[END_REF], extends our results by considering more realistic random schedules for the time instants of drug intakes: indeed, the use a homogeneous Poisson law is mathematically handy, but somewhat restrictive: in real situations, the schedule is probably much more complex. [START_REF] Fermin | Modeling patient poor compliance in the multi-IV administration case with Piecewise Deterministic Markov Models[END_REF] uses the powerful mathematical framework known as Piecewise Deterministic Markov Model to deal with general drug intake schedules, at the expense of a more complex theoretical apparatus.

We do not address in this work the problem of estimating the parameters of the various models, in particular the law of the random intake times and doses. A recent article dealing with this issue is [START_REF] Barrière | A Bayesian approach for the estimation of patient compliance based on the last sampling information[END_REF].

The remaining of this work is organized as follows: for the reader's convenience, we present an overview of our main findings pertaining to the various models and their practical implications in Section 2. In particular, we compare the variabilities in our poor compliance models to the ones in the cases of (a) full compliance of a single patient and (b) full compliance in a population where one takes into account variability due to differing elimination rates or clearances between individuals. This allows us to highlight "equivalent scenarios" where, for instance, we find the parameters of a non-compliant patient that will yield the same variance in concentration as in a whole population with given distribution of elimination rates, or, which is the same, of a single individual whose elimination rate is unknown and is modelled as a random variable. A reader not interested in the mathematical details may concentrate on Section 2 to get a quick summary of our work, and refer if needed to the details in the following sections. In Section 3, we set up the basic continuous time model for deterministic dose random time drug intake and study its variability (Section 3.1) and regularity (Section 3.2). In Section 4, we analyse the random doses version of this model. As before, we study its variability (Section 4.1) and regularity (Section 4.2). We present the discrete-time version of the deterministic dose random time model in Section 5. We derive the discrete time concentration in Section 5.1 and study its variability in Section 5.2. In Section 5.3, we show that the discrete-time model tends to the continuous-time one when the discretization step tends to 0. Section 5.4 describes the complex regularity behaviour of the probability distribution of the concentration in the discrete-time model. Finally, Section 6 deals with the discrete-time random-dose case: its variability is described in Section 6.1, and its limiting behaviour when the discretization step tends to 0 is studied in Section 6.2.

2 Overview: pharmacokinetic implication of non compliance

Purely deterministic model (full compliance)

To assess the impact of poor compliance, it is useful to contrast it with situations of full compliance. We derive in this section the variability in concentration for a single fully compliant patient. The next section will deal with the variability in a population of fully compliant patients.

Let us first recall the basic equations. We consider in this work the case of multiples oral administration and we suppose that kinetics of first order are involved. Our first model is a very classical one: it assumes that the patient takes orally a constant dose D at regularly spaced times t 0 , t 1 , . . . , t n , . . ., with t i = i λ , that is, the patient takes drugs every 1 λ units of times for some positive λ.

For completeness and future use, we briefly recall how the equations governing the absorption and elimination processes are obtained. The oral absorption process is concerned with the amount of drug at the absorption site remaining to be absorbed. The amount of drug at the absorption site, denoted A oral , is characterized by the absorption coefficient rate k a . At each drug intake time t i , A oral (t) increases by D. Between t i and t i+1 , the effect of the dose taken at t i decreases exponentially fast, with exponential speed k a . Formally:

d dt A oral (t) = ∑ i:t i ≤t Dδ(t -t i ) -k a A oral (t), (1)
with δ denoting the Dirac distribution. Thus, the rate of absorption of drug at time t is k a F A oral (t), where F is the absolute bioavailability; i.e. the fraction of each dose which is absorbed when the drug is given by the oral route. The elimination process describes by the irreversible loss of drug from the site of measurement. We assume that the drug is eliminated with a constant elimination rate k e . Denoting A central the amount of drug in the body, one gets that the rate of elimination of drug is k e A central .

The model is given from mass balance considerations: at any given time, the dose is accounted for by the amount of drug at absorption site plus the amount of drug in body plus the amount of eliminated drug. Thus, the sum of the rates of change of the drug in these compartments must equal zero, so that the rate of change of drug in body is equal to the rate of absorption minus the rate of elimination. This leads to:

d dt A central (t) = k a F A oral (t) -k e A central (t). (2)
The drug concentration at time t is defined by

C(t) = A central (t) V d , ( 3
)
where V d is the apparent volume of distribution of the drug with respect to its concentration in plasma. The solution of equation ( 1) is

A oral (t) = D ∑ i e -ka(t-t i ) 1l (t≥t i ) , (4)
where 1l (t≥t i ) denotes the indicator function of the set {t ∈ R : t ≥ t i )}, which equals 1 if t ≥ t i and 0 otherwise.

Applying the parameter-variation method to equation (2) one obtains, from equations (3) and ( 4), the following expression for the deterministic drug concentration process C d :

C d (t) = F D V d k a k a -k e ∑ i ( e -ke(t-t i ) -e -ka(t-t i ) ) 1l (t≥t i ) . ( 5
)
We will wish to compare the values of the expectation and the variance in various stochastic models to the one in the case of full compliance. Of course, in the frame of full compliance, there is no randomness involved, and one cannot define proper mean and variance. However, since the concentration varies in time, it makes sense to define the mean as the average of concentration over all values of t, and to define the variance correspondingly. In other words, we denote:

E d = lim T →∞ 1 T ∫ T 0 C d (t)dt, (6)
V ar d = lim

T →∞ 1 T ∫ T 0 (C d (t) -E d ) 2 dt (7)
for the mean and variance with respect to time of C d (t). Note that E d is closely related to the usual PK metric AU C. V ar d represent the time-average square deviation from the long term average E d . Simple computations lead to:

E d = µ e F D V d , ( 8
)
V ar d = µ e 2 ( F D V d ) 2 [ 1 1 + r + G(r) ] . ( 9
)
where

G(r) = 2 1 -r 2 ( 1 e 1/µe -1 - r e 1/rµe -1 ) -2µ e , ( 10
)
µ e := λ ke , µ a := λ ka and r := ke ka = µa µe . These quantities are related to ones considered in [START_REF] Li | A probabilistic approach for the evaluation of pharmacological effect induced by patient irregular drug intake[END_REF].

For comparison with the random models below, we note the following facts:

• When µ a → 0, we are in the case of instantaneous absorption, which is equivalent to the multi-IV case. In this frame, we recover the results of [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF].

• In Formula ( 9), the term G(r), given by Formula (10), is always negative, as may easily be checked.

• The term G(r), seen as a function of r, is increasing. It is thus minimum when r equals 0, and it tends to 0 when r tends to infinity.

• When µ e tends to 0, the variance tends to 0 at rate µ 2 e when µ a ̸ = 0 and at rate µ e when µ a = 0.

• For a fixed mean E d = 1, and when µ e tends to infinity, V ar d tends to 0 at speed 1

µ 2 e .
When µ e tend to 0, the variance tends to 1 2µa when µ a ̸ = 0, and to infinity at speed

1 µe when µ a = 0. • P(|C d -E(C d )| ≥ γE(C d )) ≤ 1
µ * 2 γ 2 for γ large enough, where µ * = max{µ e , µ a } (see Formula (25) for a more precise expression). Assume µ a , µ e tend to infinity at the same rate. Then, the bound is 1 µ * 4 γ 2 . These formulas quantify the obvious fact that the variability of concentration is a decreasing function of the number of takes per unit time. As we show below, non compliance amplifies the variability of the concentration in a way we will precisely measure.

Full compliance with population parameter variability

In the previous section, we have characterized the variability for a deterministic PK model defined by a unique set of PK parameters, i.e. the case of a single compliant subject with perfectly known pharmacokinetic parameters. We now investigate variability in a Pop-PK model for the case of a whole population. A Pop-PK involve randomly distributed PK parameter to reflect the peculiarity of each patient and other unexplained variability. A typical assumption is to consider that individual parameters are log-normally distributed in the population. The log-normal distribution is a common modelling for positive continuous quantities, and seems to fit well with the finding of large scale Pop-PK studies [START_REF] Gaudreault | A population pharmacokinetic model for the complex systemic absorption of ropivacaine after femoral nerve block in patients undergoing knee surgery[END_REF][START_REF] Sheiner | Estimation of population characteristics of pharmacokinetic parameters from routine clinical data[END_REF]. We also assume that each subject in the population is fully compliant and that their drug intake frequency λ is the same.

Let us first consider the case where the elimination rate is random and the other parameters such as D, V d , F , k a are constants. This amounts to replacing k e by k e U ke in the preceding computations, where k e is the typical parameter value of the population and U ke is a dimensionless random variable that represent the variability among subject, which is assumed to follows a log-normal distribution with parameters m ke = 0 and σ 2 ke . The choice m ke = 0 ensures that the reference value k e is the median of the distribution. We wish to examine the impact of this variability on the long-term average concentration E d , given by equations ( 6) and [START_REF] Gaudreault | A population pharmacokinetic model for the complex systemic absorption of ropivacaine after femoral nerve block in patients undergoing knee surgery[END_REF]. The variability of E d is solely due to the variation of k e U ke . In other words, we consider the variable E d (k e U Ke ) with U ke following a log-normal distribution, and measure its variability. Thus, the mean and variance of E d across the population are given by:

E pop (E d ) = λ k e F D V d E pop ( 1 U ke ) = λ k e F D V d e σ 2 ke /2 , Var pop (E d ) = ( λ k e F D V d ) 2 Var pop ( 1 U ke ) = ( λ k e F D V d ) 2 e σ 2 ke ( e σ 2 ke -1
) .

Let us rewrite these formulas in terms µ e and the coefficient of variation of the elimi-nation rate CV ke = √ e σ 2 ke -1:

E pop (E d ) = µ e F D V d √ 1 + CV 2 ke , Var pop (E d ) = ( µ e F D V d ) 2 CV 2 ke ( 1 + CV 2 ke
) .

The same approach applies if one wishes to take into account in addition the variability of V d and F across the population. In this case, k a is assumed constant and one simultaneously replaces k e by k e U ke , V d by V d U V d and F by F U F , where U ke , U V d and U F are random variables log-normally distributed with median parameters m ke = 0, m V d = 0, m F = 0 and variability parameters σ 2 ke , σ 2 V d , σ 2 F respectively. To perform the computations, the knowledge of the joint distribution of U ke and U V d , or at least that of U ke U V d is needed. Joint distributions are not usually reported in Pop-PK studies, but since the product k e V d is nothing but the clearance Cl, we can see that in this case we have:

Cl = k e U ke V d U V d = Cl 0 U Cl ,
where Cl 0 = k e V d and U Cl = U ke U V d represent the variability of clearance, whose distribution is more frequently reported in Pop-PK studies. Note that the log-normality of clearance typically reported in the literature is compatible with the joint log-normality of (U ke , U V d ). Thus, if we assume that U Cl and U F are log-normally distributed random variables with parameters m Cl = 0, m F = 0 and variability parameters σ 2 Cl , σ 2 F respectively, and we assume that the variables U Cl and U F are correlated with correlation ρ = cov(log(U ke ), log(U V d ))/σ Cl σ F ∈ [-1, 1], computations similar to the ones above lead to:

E pop (E d ) = µ e F D V d √ 1 + CV 2 , ( 11
)
Var pop (E d ) = ( µ e F D V d ) 2 CV 2 ( 1 + CV 2 ) , ( 12
)
where

CV = √ e σ 2
Cl +σ 2 F +2ρσ Cl σ F -1 is the coefficient of total variation. Note that individual variations of F are difficult to distinguish from ones of Cl. We have therefore chosen to control the variability through the coefficient of total variation CV , which encompasses the variability of k e , V d and F .

For consistency, we will always use in the sequel the following PK parameters similar to those of imatinib, see [START_REF] Barrière | A Bayesian approach for the estimation of patient compliance based on the last sampling information[END_REF][START_REF] Barrière | Compliance spectrum as a drug fingerprint of drug intake and drug disposition[END_REF][START_REF] Widmer | Population pharmacokinetics of imatinib and the role of α1-acid glycoprotein[END_REF]. They are given in Table 1 Suppose that the correlation parameter ρ is 0. Then the coefficient of total variation is CV = 44.74%, which correspond to a standard deviation of 0.4272. Assuming that doses are taken every Cl . Note that the coefficient of total variation CV is an increasing function of ρ, that reaches its maximum CV max = 65.22% when ρ = 1 and its minimum CV min = 10.31% when ρ = -1. We display on Figure 1 eight sample evolutions of the concentration for the first 14 days in this scenario. The simulations show that a steady state is quickly reached. In addition, all but one samples lie within one standard deviation

σ pop = √ Var pop (E d ) of the population mean µ pop = E pop (E d ).
In the following sections, we return to the case of a single subject with perfectly known PK parameters, but with poorly compliant behaviour, and we will compare the variability in the various cases.

Continuous model with random intake instants and deterministic doses

In the context of poor compliance the t i are not fixed, but are rather modelled as random variables. We denote these stochastic time instants as (T i ) ∈N and we assume that the intervals S i = T i -T i-1 between two doses are i.i.d. with exponential distribution of parameter λ. In other words, the sequence (T i ) i∈N is supposed to be a homogeneous Poisson process, and the mean duration between two drug intakes is equal to 1 λ . The stochastic concentration at time t reads :

C(t) = F D V d k a k a -k e ∑ i ( e -ke(t-T i ) -e -ka(t-T i ) ) 1l (t≥T i ) . ( 13
)
We illustrate this model on Figure 2, where we display eight sample paths of the evolution of C(t), simulated using the parameters given in Table 1 andλ We will show in Section 3 that C(t) has a well defined limit when t tends to infinity. We are interested in the variability and regularity properties of this steady state, or asymptotic concentration, which is denoted C. The following equations are established in Section 3.1:

-1 = 12h.
E(C) = µ e F D V d , ( 14
)
Var(C) = µ e 2 ( F D V d ) 2 1 1 + r . ( 15
)
In view of ( 14) and ( 15), we remark that:

• As expected, the means are the same in the deterministic and random models.

• As is intuitively obvious, the variance in the case of full compliance is always smaller than the one in the random situation. This can be observed from the fact that the term

G(r) = 2 1-r 2 ( 1 e 1/µe -1 - r e 1/rµe -1 )
-2µ e in Formula ( 9) is always negative. Formulas ( 9) and ( 15) quantify in a precise way how much variability is increased as an effect of poor compliance.

• When µ e tends to 0, the variance tend to 0, with the same speed as in the deterministic model.

• For a fixed mean equal to 1, when µ e tends to infinity, the variance tends to 0 at speed 1 µe . This is to be compared to the faster speed of convergence of 1 µ 2 e in the deterministic case. When µ e tend to 0, the variance tends to 1 2µa for µ a ̸ = 0, which is strictly larger than the corresponding limit in the deterministic case. Finally, when µ a = 0, the variance tends to infinity at same rate as in the deterministic model ( 1µe ).

• P(|C -E(C)| ≥ γE(C)) ≤ 1
2µ * γ 2 for γ enough large and for all µ * (recall that µ * = max{µ e , µ a }). This is a slower speed of convergence than in the deterministic case, which is 1 µ * 2 . • The probability that the long-term concentration exceeds a given level γ decays as

1 γ 2 . More precisely, P(C ≥ γ) ≤ 7 6γ 2 ( F D V d ) 2 µ 2 e ( 1 + 1 2(µe-µa)
)

.

If E(C) = µ e F D V d
is assumed to be constant, one sees that this probability decays at the same speed for all values of µ e , µ a .

The coefficient of variation of the limit concentration C is 1/ √ 2(µ e + µ a ); as a consequence, variability is larger for smaller values of µ * , i.e. when the absorption process or the elimination process is fast compared to the mean frequency of drug intake. In Section 3.2, we will show that µ * also controls the regularity of the limit concentration distribution.

The difference between the variance V ar d and Var(C), given by equations ( 9) and ( 15) respectively, is largest when r is small, and vanishes when r tends to infinity. The interpretation of this fact is straightforward: for a fixed k e , a large r means a small k a and thus a slow absorption. In this case, the effect of a drug intake on the system takes a very long time to appear, and thus a random delay in the intake has almost no effect. On the contrary, a small r translates into fast absorption, and any irregularity in the schedule has a noticeable impact. We also remark that, because of the "damping" introduced by k a in the oral model, the increase of variance due to random intakes is smaller here than in the intravenous model studied in [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF]. When r tends to 0, we recover the results of [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF] since the damping vanishes. Again, the formulas above quantify precisely these effects.

In a way similar to what was observed in the deterministic model, the variability of concentration is a decreasing function of the expected number of takes per unit time λ -1 : as is intuitively clear, increasing the mean frequency of intakes while keeping constant the average quantity of administrated drug diminishes the negative impact of poor compliance in terms of the probability of departing significantly from the mean concentration. In order to compare the variances in the deterministic [START_REF] Kingman | Poisson Processes[END_REF] and non compliant cases [START_REF] Pratsiovytyi | A class of random variables of the Jessen-Wintner type[END_REF], we plot their evolutions as a function of µ e and µ a in two situations. First, in Figure 3, we let µ e and µ a vary and keep the other parameters, i.e. D, F and V d, constants. This corresponds for instance to the case where the number of doses per unit time (or average number of doses per unit time in the stochastic case) evolves, while maintaining everything else unchanged. As can be seen on Figure 3, when µ e tends to infinity, the variance in the deterministic case reaches a plateau, while the random situation leads to an unbounded variance. Hence non-compliance has here a dramatic effect. When µ e tends to 0, both variances tend to 0, at rate µ 2 e when µ a ̸ = 0 and at rate µ e when µ a = 0. The conclusions in this frame are however somewhat unrealistic, since they lead to an unbounded increase of the mean when µ e tends to infinity. This is why, in Figure 4, we consider the more reasonable situation where µ e and µ a vary but the mean is kept constant. From equations ( 9) and [START_REF] Pratsiovytyi | A class of random variables of the Jessen-Wintner type[END_REF], one sees that this simply translates into ensuring that µ e F D V d is constant. Thus we are in the case where, for instance, one increases the frequency of drug intake while decreasing accordingly the unit dose. Here, we take µ e F D V d = 0.8743, which corresponds to the parameters given in Table 1, and λ -1 = 12h. Figure 4 indicates that both variances tend to 0 when µ e tends to infinity. The speed of convergence is however faster, as expected, in the deterministic case ( 1

µ 2 e
) than in the random one ( 1 µe ). Consider now the case where µ e tend to 0: when µ a ̸ = 0, the variance in the random frame tends to 1 2µa and to a strictly smaller value in the deterministic case. When µ a = 0, both variances tend to infinity at same rate ( 1 µe ). For instance, the point (1/24; 0.079) in the graph on Figure 5 means that a compliant patient taking a dose every day will have same concentration variability as a non compliant one taking in average a dose every 12.66 = 1/0.079 h, when the mean concentration, mean dose per day and all other parameters are the same in both situations.

Let us now compare the variabilities induced by non compliance in the random instants model and by differing PK-parameters in a fully compliant population as studied in Section 2.2. We consider the case where k a and V d are constant and k e , F are replaced by k e U ke and F U F , with U ke and U F log-normally distributed random variables with median parameters m ke = 0, m F = 0 and variability parameters σ 2 ke , σ 2 F respectively. As in Figure 5, we plot on Figure 6 the values of λ in the random model as a function of σ 2 ke and σ 2

F that yield the same variance of concentration when the mean of concentration, the mean dose per day and the other parameters are the same in both situations. From equations ( 11), ( 12), ( 14) and ( 15), one sees that this amounts to setting:

λ = k e 2(1 + r) ( e σ 2 ke +σ 2 F -1
)

,
where λ is the mean number of doses per hours of the non compliant patient, k e is both his elimination rate and the mean elimination rate in the population, σ 2 ke and σ 2 F are the variances of the elimination rate and the bioavailability in the population respectively. For example, the point (0.12, 0.06, 0.1) in the graph on Figure 6 has the following meaning: the variance of the concentration for a population of compliant patients with σ 2 ke = 0.12 and σ 2 F = 0.06 is the same as the one of a single non compliant patient taking in average a dose every 10 (= 1/0.1) hours. Likewise, a single non compliant patient taking in average a dose every day displays the same variance in concentration as a population of fully compliant individuals with σ 2 ke + σ 2 F = 0.38, i.e. a coefficient of total variation equal to 68%. Note that the same conclusions hold if one consider a population where the clearance instead of the elimination rate follows a lognormal distribution.

Let us now consider the regularity of the distribution of C. It is shown in Section 3.2 that the smoothness of the distribution of the concentration is governed by µ * = max{µ e , µ a }. We prove that, in the long term, the cumulative probability distribution of drug concentration may display two types of behaviours: when µ * is larger than one, the distribution is regular, while, when it is smaller than one, it is singular only at the origin. Thus, µ * = 1 is the critical value below which the moments of intakes are too scarce with respect to the elimination and absorption rates, resulting in a high probability of having too small a concentration of drugs.

To illustrate this fact, we plot on Figure 7 histograms representing the empirical probability distribution of C(T ), for a fixed time T , in two particular cases: µ * = 2 and µ * = 0.5. The value of D in each case was adjusted accordingly in order to keep a constant mean. These histograms were obtained by simulating 50, 000 independent sample paths of concentration in each scenario until time T and distributing the outcomes into 100 evenly spaced bins. The time T was chosen large enough so that the steady state has been reached (we have set T = 100 max{1/λ, 1/k e }). The singularity of the distribution of C at the origin in the case µ * < 1 (second histogram) manifests itself through the sharp spike in the first bin. The reader is referred to Section 3 for more details on this model.

Continuous model with random intake instants and random doses

In this section, we sum up the results obtained in Section 4, where we generalize the continuous model of Section 2.3 to allow for random doses. The idea is that the careless patient that has an irregular schedule is likely to also mess with the doses. For instance, he might take a double dose to make up for a missing one. Formally, this translates into the fact that the quantity D in ( 5) is now a random variable that may vary at each take rather than being a constant. In other words, the concentration, denoted C rd (t), is now given by:

C rd (t) = F D V d k a k a -k e ∑ i D i ( e -ke(t-T i ) -e -ka(t-T i ) ) 1l (t≥T i ) , ( 16
)
where D, F and V d are constants, the (T i ) i∈N again form a homogeneous Poisson process, and the (D i ) i∈N are random variables. The dose taken at time T i is DD i . It thus seems natural to assume that E(D i ) = 1 (i.e on average, the patient takes the required dose), that D i is supported on R + , and that it has compact support, i.e. the patient cannot take arbitrarily large doses (although we shall need only weaker assumptions). We illustrate this model by showing on Figures 8 and9 simulated sample paths of the concentration in two particular cases considered in more details in Section 4: in the first one, the random factors D i follow a uniform distribution on an interval [a, b], while in the second one they follow a discrete distribution taking two possible values d 1 and d 2 with probability q 1 , q 2 . The PK parameters chosen for the simulation of the sample paths are the ones given in Table 1 with an average time between intakes equal to λ -1 = 12h. The parameters for the random doses in both cases are given in Table 2. We prove in Section 4.1 that the mean value in this model is the same as in the previous ones, i.e.

d 1 = 0.4 d 2 = 1.9 q 1 = 0.6 q 2 = 0.4 a = 0.2 b = 1.8
E(C rd (t)) = E(C(t)).

Furthermore, under the assumption that the doses (D i ) i∈N do no depend on the times (T i ) i∈N , we show the following equality:

Var(C rd (t)) = (1 + Var(D 1 ))Var(C(t)).
In other words, the variance in the case of random dosing is simply the one of the deterministic dose case multiplied by 1 + Var(D 1 ). As a consequence, estimating quantities such as P(|C -E(C)| ≥ γE(C)) or P(C ≥ γ) or comparing to the case of (population) full compliance is readily done once Var(D 1 ) is known. In Section 4.1, we detail two particular cases: the first one is where D i is uniformly distributed on a interval, and the second one is where it takes values in a discrete set. We characterize the situations leading to the largest variances for both distributions (formulas (33) and ( 35)).

In contrast to variability, random dosing does not affect the regularity of the distribution of the long term concentration: this distribution is again smooth when µ * > 1 and is singular at 0 otherwise.

Discrete model with random intake instants and deterministic doses

We study a time-discretized version of the random model presented in Section 2.3. The general idea is that, instead of taking the drug at arbitrary time instants t ∈ R + , the patient will only do so at times which form a random subset of {t n l = l/n, l ∈ N}, where n is a fixed number. There are two main reasons for considering such a model. First, there are indeed natural situations where a time discretization does occur. For instance, the medication must sometimes be taken at precise moments, like before lunch. Many people will always have their lunch at a fixed time, like certain workers or many older people. Second, as we show in Section 5.3, the time discretized model tends to the continuous one when the discretization step tends to 0. Thus, for n large enough, the practical difference between both models vanishes. Nonetheless, the discrete model displays various interesting and intriguing features that are not present in the continuous one. Let us finally mention that considering the discretized model is very close to sampling the concentration of the continuous one. Since blood concentrations cannot be monitored continuously, the outcome of any clinical study is discrete in nature, which gives further justification for the discrete model.

We show in Section 5.1 that the steady state discretized concentration Y (p) in the situation where the patient takes a dose D with probability p independently at each time j/n reads:

Y (p) = F D V d k a k a -k e ∞ ∑ j=0 ( β j e -β a j ) X j , ( 17
)
where the (X j ) j∈N are i.i.d. Bernoulli random variables with parameter p (i.e. X j = 0 with probability 1 -p and X j = 1 with probability p). The parameter β e = e -ke/n is the elimination rate for one time step while β a = e -ka/n is the absorption rate for one time step. As explained in Section 5.1, in the discrete model, the number of drug intakes per unit time is np, thus one has to set p = λ n to ensure correspondence with the continuous model. Set α := F D V d ka ka-ke . The mean and variance of the discretized concentration are given by the following formulas:

E disc := E[Y (p) ] = αp ( 1 1 -β e - 1 1 -β a ) , V ar disc := Var(Y (p) ) = α 2 p(1 -p) ( 1 1 -β 2 e - 2 1 -β e β a + 1 1 -β 2 a )
.

Note that E disc and V ar disc tend respectively to the mean and variance in the continuous model with random instants of intakes when p tends to 0.

If ones fixes E disc = 1, then the variance in the discrete model becomes:

V ar disc = 1 -p p (1 -β e )(1 -β a )(1 + β e β a ) (1 + β e )(1 + β a )(1 -β e β a )
.

Thus V ar disc tends to 0 at speed 1-p 2µe when µ e tends to infinity, which is the same rate as V ar, the variance in the continuous random model of Section 3.1. When µ e tends to 0, V ar disc tends to 1-p p 1-βa 1+βa in contrast to both the deterministic and random continuous cases, where the variance tends to infinity. Figure 10 displays the evolution of the ratio V ar disc /V ar as a function of µ e and µ a when the mean is kept constant.

Formula (36) in Section 5.2 gives the probability that the concentration departs significantly from its mean. It is interesting to compare it with (24) to contrast the impact of non-compliance in the discrete and continuous time models, and to (25) to measure the added variability with respect to the fully compliant situation.

The study of the regularity of the distribution is much more involved in the discrete case than in the continuous one, as the former may display a fractal behaviour. Indeed, for some values of the parameters, the distribution will be everywhere singular. This means that the evolution of the distribution varies wildly for most values of the concentration. In addition, the dichotomy smooth/singular is not governed by µ * alone, but by complex relations between k e , k a and the discretization step. More precisely, the relevant parameter here is β * = β e ∨ β a and the regularity of the discretized model is very different depending on whether

β * < 1 2 or β * ≥ 1 2 . When β * < 1
2 , the distribution of Y (p) is singular with respect to the Lebesgue measure. This means that, in this situation, the probability distribution of the concentration will be highly irregular : it will vary erratically, taking only very particular values, and the probability that it ranges in some interval varies wildly when the interval changes. Note that β * < 1 2 is equivalent to p > µ * log(2) or k * > n log(2) with k * = k e ∧ k a . This reflects the fact that, for fixed k * , the time instants at which the patient is supposed to take his drugs are sufficiently spaced, and that he forgets to do so randomly at some of these instants.

The case β * > 1 2 is mathematically much more delicate and is not completely settled. What is known is that the distribution of Y (p) is absolutely continuous for almost all β * in ( 1 2 , 1) when p = 1 2 : in effect, this means that the distribution of the concentration has a probability density and thus we are in the usual situation where the probability that the concentration ranges in some interval varies smoothly when the interval changes.

See Section 5.4 for more details.

Discrete model with random intake instants and random doses

Similarly to Section 2.4, we consider a random dose variant of the previous model. The steady state concentration now reads:

Y rd (p) = α ∞ ∑ j=0 (β j e -β j a )D j X j , ( 18
)
which accounts for the fact that the stochastic dose at time j/n is equal to 0 or DD j , where the random variables D j are i.i.d. with mean 1 and take values in an interval

[d m , d M ] with 0 < d m < d M < ∞.
The mean E rd disc and variance V ar rd disc in the random-dose discretized model are given by (see Section 6):

E rd disc = αp ( 1 1 -β e - 1 1 -β a ) = E disc , V ar rd disc = ( Var(D 1 ) 1 -p + 1 ) V ar disc .
Random-dosing thus results in multiplying the variance of the deterministic-dose discretized model by Var(D 1 )

1-p + 1. It is then easy to study the variability of the random-dose discretized model using the analysis developed in Section 5.2 for the discretized model.

We do not have at this time any significant result concerning the regularity of the distribution of the concentration in this model. This ends the overview of our models and their main features. The reader will find in the next sections the precise derivations of the results mentioned above, along with a more detailed analysis of the models.

Continuous model with random intake instants and deterministic doses

In this section we study the main properties of the random concentration given by ( 13), which we recall here for convenience:

C(t) = F D V d k a k a -k e ∑ i ( e -ke(t-T i ) -e -ka(t-T i ) ) 1l (t≥T i ) ,
where (T i ) i∈N is a homogeneous Poisson process with parameter λ. The probability distribution of C(t) may be described through its characteristic function φ t . It is easily obtained by applying Campbell's theorem [START_REF] Kingman | Poisson Processes[END_REF]. In our case, this yields:

φ t (θ) = exp { λ ∫ t 0 ( exp { iθα ( e -ke(t-x) -e -ka(t-x) )} -1 ) dx } ×exp { iθα ( e -ket -e -kat
)} .

where we recall that α := F D V d ka ka-ke . The change of variable u = e -(t-x) leads to:

Proposition 1. The characteristic function of C(t) is φ t (θ) = exp { λ ∫ 1 e -t exp { iθα ( u ke -u ka )} -1 u du + iθα ( e -ket -e -kat ) } . (19)

Variability of the concentration

The characteristic function allows us to compute the mean and variance of C(t):

E(C(t)) = λα ( 1 -e -ket k e - 1 -e -kat k a ) + α ( e -ket -e -kat
) , so that, in the long term,

lim t→∞ E[C(t)] = µ e F D V d , ( 20
)
Likewise, the variance of C(t) reads:

Var(C(t)) = λα 2 ( 1 -e -2ket 2k e + 1 -e -2kat 2k a -2 1 -e -(ka+ke)t k a + k e ) ,
from which one gets

V ar = lim t→∞ Var(C(t)) = µ e 2 1 1 + r ( F D V d ) 2 , ( 21
)
with r := ke ka . Note that the convergence in ( 20) and ( 21) are exponential: only a few cycles are needed before the steady state is reached. The same remark applies to all convergences below.

In the next section, we show that, when t tends to infinity, C(t) also converges to a well defined random variable, denoted C, and we investigate in details some of its properties. Before, let us give a final result of interest pertaining to the variability of the concentration in the non-compliant case. More precisely, the following proposition yields bounds on the probability that the concentration exceeds a given (large) level, or departs significantly from its mean. Proposition 2. For γ large enough,

P(C(t) ≥ γ) ≤ 7 6γ 2 ( F D V d ) 2 ( µ e µ e -µ a ) 2 H(t, k e , k a , µ e , µ a ), (22)
P(C ≥ γ) ≤ 7 6γ 2 ( F D V d ) 2 µ 2 e ( 1 + 1 2(µ e -µ a )
)

. ( 23)

P(|C -E[C]| ≥ γE[C]) ≤ 1 2µ e γ 2 1 1 + r . ( 24
)
where

H(t, k e , k a , µ e , µ a ) = [ (µ e -1) ( 1 -e -ket ) -(µ a -1) ( 1 -e -kat )] 2 + [ µe
Proof. This is a direct application of the classical bound (see, e.g. [START_REF] Loeve | Probability Theory I, Fourth Edition[END_REF], p. 209):

P(C(t) ≥ γ) ≤ 7γ ∫ 1 γ 0 (1 -Re(φ t (θ))) dθ,
valid for γ > 0, and where Re denotes the real part. Indeed,

Re(φ t (θ)) = cos ( θα(e -ket -e -kat ) + λ ∫ 1 e -t sin(iθα(u ke -u ka )) u du ) × exp { λ ∫ 1 e -t cos(θα(u ke -u ka ))-1 u du } ,
and routine estimates yield (22). Inequality (23) follows in a similar way. Finally, ( 24) is simply Chebychev inequality.

Note that, in the deterministic case (full compliance), and with the same definition of the variance given by ( 7) in Section 2.1, one has in place of (24):

P(|C d (t) -E[C d ]| ≥ γE[C d ]) ≤ 1 2µ e γ 2 [ 1 1 + r + G(r) ] , (25)
with G(r) given by [START_REF] Iskedjian | Relationship between daily dose frequency and adherence to antihypertensive pharmacotherapy: evidence from a meta-analysis[END_REF]. This is another quantitative way to measure by how much the probability of differing from the mean will be larger in the non-compliant case. For instance, when µ e tends to infinity but µ a remains bounded (which implies that r tends to 0), the above bound is of the order of 1 µ 2 e γ 2 , and thus much smaller than the one in (24). When µ a tends to infinity but µ e remains bounded (i.e. r tends to infinity), the bound is of the order of 1 µ 2 a γ 2 , again much smaller than the one in (24). Finally, when µ e and µ a tend to infinity at the same rate, the bound in the deterministic case is of the order of 1 µ 4 e γ 2 . This was illustrated on Figure 3.

Regularity of the limit distribution

In this section, we study the long term behavior of the drug concentration, that is, the distribution function of the limit C := lim t→∞ C(t), where the limit is taken in the sense of convergence in distribution.

Proposition 3. The random variable C(t) converge in distribution, when t tends to infinity, to a well define random variable C whose characteristic function is

φ(θ) = exp { λ ∫ 1 0 exp { iθα ( u ke -u ka )} -1 u du } . ( 26
)
Proof. When t tends to infinity, φ t tends pointwise to φ, which is continuous at θ = 0. By Lévy's theorem, this implies the result.

22

Note that the distribution of C invariant by time reversal: looking "backwards" in time, one see that C(t) has the same law as

C ′ (t) = F D V d k a k a -k e ∑ i ( e -keT i -e -kaT i ) 1l (t≥T i ) ;
as t tend to infinity, the random variables C ′ (t) converge almost surely to

C ′ = F D V d k a k a -k e ∑ i ( e -keT i -e -kaT i
) ,

which therefore has the same distribution as C. In the sequel, we shall write

C = F D V d k a k a -k e ∑ i ( e -keT i -e -kaT i
) ,

since we are only interested in distributional properties.

Proposition 4. The characteristic function φ satisfies

|φ(θ)| ∼ K|θ| -µ * , when θ → ∞,
with K a positive constant and µ * = λ (ke∧ka) = max{µ e , µ a }. Proof. One computes:

|φ(θ)| = exp { -λ ∫ 1 0 1 -cos ( θα ( u ke -u ka )) u du } =: exp {-λI(θ)} .
Set 0 < δ < 1, and decompose the integral I(θ) as follows:

I(θ) = I 1 (θ) -I 2 (θ) -log(δ), (27)
where

I 1 (θ) = ∫ δ 0 1 -cos ( θα ( u ke -u ka )) u du, I 2 (θ) = ∫ 1 δ cos ( θα ( u ke -u ka )) u du. ( 28
)
When 0 ≤ u ≤ δ, we have that α ( u ke -u ka ) ∼ |α|u ke∧ka , and thus

I 1 (θ) ∼ ∫ δ 0 1 -cos ( θαu ke∧ka ) u du ∼ 1 k e ∧ k a log(θ|α|) + γ e k e ∧ k a + log(δ) + O ( 1 θ ) , (29)
where γ e is the Euler constant.

Denote

h(u) = α ( u ke -u ka )
. The function h has a global maximum at u 0 = ( ke ka

)1

ka-ke

with h " (u 0 ) < 0. Stationary phase arguments imply that, when θ tends to infinity,

I 2 (θ) ∼ Re ( ∫ 1 δ e iθα(u ke -u ka ) u du ) ∼ Re ( e i(θh(u 0 )-π/4) √ 2π u 0 √ θ|h " (u 0 )| ) = cos(θh(u 0 ) -π/4) u 0 √ 2π θ|h " (u 0 )| = O ( 1 √ θ
) .

(30) Formulas ( 27), ( 28), ( 29) and (30) entail that

I(θ) = 1 k e ∧ k a log(θ|α|) + γ e k e ∧ k a + O ( 1 √ θ
) .

This implies that there exists K > 0 such that |φ(θ)| ∼ K|θ| -µ * when θ tend to infinity.

We denote by F the probability distribution function of C, associated to φ. Then, from Proposition 4 we have the following results with respect to the regularity of F :

1. F has an L 2 density 1 if and only if µ * > 1 2 . 2. For µ * < 1 and 0 < ε < 1, ∫ ∞ 1 θ µ * -ε |φ(θ)|dθ < ∞ thus F ∈ Lip(µ * -ε). 3. For µ * < 1 2 , 1 T ∫ T -T |φ(θ)| 2 dθ = O ( T -2µ * ) thus F ∈ Lip(µ * ).
4. A classical Tauberian theorem (see, e.g. [START_REF] Feller | An introduction to probability theory and its applications[END_REF], p.445) entails that:

F (ε) ∼ e -µ * γe |α| -µ * Γ(µ * + 1) |ε| -µ * ,
when ε → 0. Then, F is not differentiable at 0 when µ * < 1 and it has a finite non vanishing derivative at 0 exactly when µ * = 1. [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF], we have that, for any x > 0,

From Proposition 3 in

F (x + ε) -F (x) = O(ε),
when ε → 0 + . This implies that 0 is the only possibly singular point of F .

The practical meaning of these results is that, when λ < min{k a , k e } the instant of intakes are too scarce, with respect to the absorption and the elimination rates, resulting in a high probability of having too small concentration of drugs.

Continuous model with random intake instants and random doses

In this section, we consider the generalization of the continuous model of Section 3 to allow for random doses. We study the main properties of the random concentration C rd (t), given by ( 16), which we recall here for convenience:

C rd (t) = F D V d k a k a -k e ∑ i D i ( e -ke(t-T i ) -e -ka(t-T i ) ) 1l (t≥T i ) ,
where the (T i ) i∈N again form a Poisson process, and the (D i ) i∈N are random variables. F and V d are constants. At time T i , the dose taken is DD i , and we assume that E(D i ) = 1, that D i is supported on R + , and that it has compact support (although we shall need only weaker assumptions).

The process C rd (t) thus defined is a marked Poisson process. In this work, we shall assume that the (D i ) i∈N are independent and identically distributed random variables, where each D i may depend on T i but is independent from the (T j ) j̸ =i . This makes sense from a pharmacokinetic point of view, since it seems plausible that the patient will not adjust his dose at time T i on the basis of his past or future behavior except for the time lag from the previous take, although it would maybe be desirable to let D i depend also on D i-1 . We denote by ν(T, .) the conditional distribution of D i knowing that T i = T . Our assumptions allow to apply a generalized form of Campbell theorem [START_REF] Kingman | Poisson Processes[END_REF] to the effect that the characteristic function φ rd t of C rd (t) is given by:

φ rd t (θ) = exp { λ ∫ t 0 ∫ A ( e iθαuh(t-x) -1 ) ν(x, du)dx + iθαh(t) } (31)
where, as before, α = F D V d ka ka-ke , h(t) = e -ket -e -kat and A is the support of the D i .

Variability of the concentration

From (31), one deduces easily the mean and variance of C rd :

(φ rd t ) ′ (θ) = φ rd t (θ) [ λ ∫ t 0 ∫ A iαuh(t -x)e iθαuh(t-x) ν(x, du)dx + iαh(t) ] ,
and thus

(φ rd t ) ′ (0) = iα ( µ e (1 -e -ket ) -µ a (1 -e -kat ) ) E x (D 1 ) + iα ( e -ket -e -kat
) .

where E x (D 1 ) = ∫ A uν(x, du) is the expectation of D i knowing that T i = x, which is equal to one for any x by assumption. Thus we find that, not surprisingly:

E(C rd (t)) = E(C(t)).
Likewise,

(φ rd t )"(θ) = φ rd t (θ) [ ( λ ∫ t 0 ∫ A iαuh(t -x)e iθαuh(t-x) ν(x, du)dx -α 2 h 2 (t -x) ) 2 -λ ∫ t 0 ∫ A α 2 u 2 h 2 (t -x)e iθαuh(t-x) ν(x, du)dx ] .
Using that, by definition,

∫ A u 2 ν(x, du) = E x [D 2 1 ]
, this entails:

(φ rd t ) ′′ (0) = -E[C(t)] 2 -λα 2 ∫ t 0 h 2 (t -x) ∫ A u 2 ν(x, du)dx = -E[C(t)] 2 -λα 2 ∫ t 0 h 2 (t -x)E x (D 2 1 )dx or, since (φ rd t ) ′′ (0) = -E[C(t)] 2 -Var(C rd (t)), Var(C rd (t)) = λα 2 ∫ t 0 h 2 (t -x)E x [D 2 1 ]dx.
However,

E x (D 2 1 ) = Var x (D 1 ) + E x [D 1 ] 2 = Var x (D 1 ) + 1 (Var x (D 1 ) denotes the variance of D i knowing that T i = x) and λα 2 ∫ t 0 h 2 (t -x)dx = Var(C(t)), thus Var(C rd (t)) = Var(C(t)) + λα 2 ∫ t 0 h 2 (t x )Var x (D 1 )dx.
Assuming that Var x (D 1 ) = Var(D 1 ), i.e. the variance of D i does not depend on the value of T i , one gets:

Var(C rd (t)) = E[D 2 1 ]Var(C(t)) = (1 + Var(D 1 ))Var(C(t)).
Thus random-dosing results in multiplying the variance of the deterministic-dose case by E(D 2 1 ). It is then easy to obtain inequalities similar to (22), ( 23) and (24) for the randomdose case.

Two particular cases may be of special interest:

• Discrete distribution: the D i assume the values {d 1 , . . . , d m } ⊂ R + with probabilities q 1 , . . . , q m , independently of the T i . We denote by C rd,d t and φ rd,d t the concentration and characteristic function. Then:

φ rd,d t (θ) = m ∏ j=1 φ t (q j λ, d j θ), ( 32
)
where we have written φ t (q j λ, •) instead of φ t (.) for the characteristic function in the deterministic-dose case, given in equation ( 19), but with parameter q j λ instead of λ. This equality also holds when t tends to infinity, yielding the long term concentration characteristic function. The variance is given by:

Var(C rd,d (t)) = E[D 2 1 ]Var(C(t)) =   m ∑ j=1 q j d 2 j   Var(C(t)).
It is of interest to characterize the situation giving the largest variance among all admissible random dosings with arbitrary m, q i , and d i . In other words, we look for the value of V max := max{Var(C rd,d (t))} subject to m > 1, (q 1 , . . . , q m ) ∈ [0, 1] m with at least two q i non zero and

∑ m j=1 q i = 1, (d 1 , . . . , d m ) ∈ [a, b] with b > a > 0, and ∑ m j=1 q i d i = 1.
It is easily shown, for instance using Lagrange multipliers, that the maximum is reached for m = 2 and d 1 = a, d 2 = b. In this case,

V max = (a + b -ab)Var(C(t)). ( 33
)
Thus the worst-case situation is when the patient "oscillates" between two dosings whose average is the prescribed one.

• Uniform distribution: the D i are uniformly distributed over [a, b] ⊂ R + , independently of the T i . We denote by C rd,u t and φ rd,u t the concentration and characteristic function. In this case, one computes:

φ rd,u t (θ) = exp { λ b -a ∫ b a ∫ 1 e -t e iθyα(u ke -u ka ) -1 u dudy + iθαh(t) } . ( 34
)
This last function is easily seen to be convergent when t tends to infinity, which gives us the long term concentration characteristic function. Note that φ rd,u tends to φ when the couple (a, b) tends to (D, D): the concentration with random doses uniformly distributed on [a, b] tends in law to the concentration with fixed dose D.

The variance is given by:

Var(C u (t)) = E(D 2 1 )Var(C(t)) = a 2 + ab + b 2 3 Var(C(t)). (35)
We can see that the choice of [a, b] maximizing the variance under the constraint E(D 1 ) = 1 is a = 0, b = 2, as expected.

Regularity of the limit distribution

We first consider the two cases of discrete and uniform random dosing.

• Discrete distribution From (32) and Proposition 4, one easily sees that:

|φ rd,d (θ)| ∼ m ∏ j=1 (d j θ) -q j µ * = Kθ -µ *
when θ tends to infinity, where K is a constant.

• Uniform distribution In this case, the limit of the function φ rd,u t , defined in Formula (34), when t tends to infinity is:

φ rd,u (θ) = exp { λ b -a ∫ b a ∫ 1 0 e iθyα(u ke -u ka ) -1 u dudy } .
The modulus of φ rd,u is

|φ rd,u (θ)| = exp { -λ (b -a) ∫ b a ∫ 1 0 1 -cos(θyα(u ke -u ka )) u dudy } , = exp { -λ (b -a) ∫ b a I(θy)dy } .
The proof of Proposition 4 shows that, for large enough θ,

I(θy) ∼ 1 k e ∧ k a log(θy|α|) + constant.
This implies that there exists K > 0 such that |φ(θ)| ∼ K|θ| -µ * when θ tends to infinity. This fact will be justified in the proof of Proposition 5.

Thus in both particular cases above, we recover the same behaviour for the characteristic function as in the deterministic-dose case. As the next proposition shows, this is in fact a general feature of all random dosings provided ν(x, du) does not depend on x, with an additional mild condition: Proposition 5. Assume that ν(x, du) = ν(du) and that ∫ A log(y)ν(dy) < ∞. Then

|φ rd (θ)| ∼ K|θ| -µ *
when θ tends to infinity, where K is a positive constant.

Proof. Thanks to the assumption on ν, one computes

|φ rd,u (θ)| = exp { -λ ∫ A ∫ 1 0 1 -cos(θyα(u ke -u ka )) u duν(dy) } = exp { -λ ∫ A I(θy)ν(dy)
} .

The proof of Proposition 4 shows that, when θ tends to infinity, -λI(θy) ∼ K -µ * (log(θ) + log(y))

for a certain constant K, with in addition -λI(θy) -K + µ * (log(θ) + log(u)) tending to 0 with a rate of convergence 1 √ θ . Thus, using the assumption on the logarithmic moment of ν, -λ

∫ A I(θy)ν(dy) ∼ K -µ * log(θ) -µ * ∫ A log(y)ν(dy)
and one finishes up the proof with the help of a dominated convergence argument to show that the difference between the right-hand side and the left-hand side in the equivalent above indeed tends to 0 when θ tends to infinity.

Proposition 5 shows that random dosing, at least when the distribution of the D i is independent of the T i , does not alter the regularity of the distribution of the long-term concentration as compared to the deterministic-dose case.

Discrete model with random intake instants and deterministic doses

We study in this section a discretization in time of the model above. In other words, instead of taking the drug at arbitrary time instants t ∈ R + , we assume that the patient will only do so at (random) times which form a subset of {t n l = l/n; l ∈ N}, where n is a fixed number. We shall first rewrite the drug concentration in this discrete setting. It will appear that the discretized concentration C d has the same law as an object that has been thoroughly studied in mathematics under the name of infinite Bernoulli convolution. We will study the variability of the discretized concentration. Then we will show that, when n tends to infinity, the discretized model indeed tends distribution-wise to the continuous-time one. We will finally study the regularity of the long term behaviour of the discretized concentration for n fixed or tending to infinity, and show that it is, under certain circumstances, singular.

The discretized concentration

The discretized model taking the following form. Let h = 1 n be the discretization step. Thus, the drug intakes can only occur at times t j = jh with j ∈ N.

From a general point of view, the discrete analog of the Poisson process is the Bernoulli process. Indeed, in the continuous framework, the Poisson process is the only counting process which has stationary and independent increments. Likewise, the only discrete counting process with the same property is the Bernoulli process. In terms of waiting times, this amounts to replacing the i.i.d., memoryless, exponential random variables S n = T n+1 -T n , by i.i.d. random variables following a geometric distribution (recall that the geometric distribution is the only memoryless discrete distribution).

We are thus led to consider a sequence (X j ) j∈N of i.i.d. Bernoulli random variables with parameter p such that X j = 1 if the patient takes the drug at time t j and X j = 0 if not.

In this discrete model the number of drug intakes per unite time is p h , so one has choose p = λh. Note that, for the model to make sense, p must be smaller than one, which translates into λh < 1.

At a fixed time t n the contribution of the j-th intake to the current concentration is

α ( β n-j e -β n-j a ) X j ,
where β e = e -keh and β a = e -kah : β e is the elimination rate for one time step and β a is the absorption rate for one time step. Thus, the total concentration at time t n is given by

C n = α n ∑ j=0 ( β n-j e -β n-j a ) X j .
Since the random variable (X j ) j∈N are independent, they are exchangeable, and in particular the vector (X 0 , . . . , X n ) has the same distribution as the reversed vector (X n , . . . , X 0 ). Hence C n is equal in distribution to

Y n = α n ∑ j=0 ( β j e -β a j ) X j .
The sequence (Y n ) converges almost surely to the random variable Y (p) , given by equation (17) in Section 2.5; i.e.,

Y (p) = α ∞ ∑ j=0 ( β j e -β a j ) X j .
The distribution of Y (p) is thus the one of the long term concentration in this model. Since Y (p) is an infinite sum of independent Bernoulli random variables, its law is an infinite convolution of Bernoulli distribution, hence its name "infinite Bernoulli convolution".

From the independence of (X j ) j∈N we have that the characteristic function of Y (p) is given by

φ p (θ) = ∞ ∏ j=0 [ (1 -p) + pe iθα(β j e -β j a )
] .

Variability of the discretized concentration

One computes easily the mean E disc and variance V ar disc in the discretized model:

E disc := E[Y (p) ] = αp ( 1 1 -β e - 1 1 -β a ) , V ar disc := Var(Y (p) ) = α 2 p(1 -p) ( 1 1 -β 2 e - 2 1 -β e β a + 1 1 -β 2 a
)

.

As we shall see in the next section, the discrete model tends in law to the continuous one when p tends to 0. For now, we simply remark that, indeed, E disc tends to µ e F D V d and V ar disc tends to µe

( F D V d ) 2
when p tends to 0, i.e. the mean and variance of (20) and

(21) (note that, when p tends to 0, β e ∼ 1 -p µe and β a ∼ 1 -p µa so β e and β a tend to 1). For E disc constant, say E disc = 1, the variance in the discrete model reads:

V ar disc = 1 -p p (1 -β e )(1 -β a )(1 + β e β a ) (1 + β e )(1 + β a )(1 -β e β a ) .
For a fixed p, we see that V ar disc tends to 0 at the speed 1-p 2µe when µ e tends to infinity, which is the same rate as V ar, the variance in the continuous random model of Section 3.1.

When µ e tends to 0, V ar disc tends to 1-p p 1-βa 1+βa in contrast to both the deterministic and random continuous cases, where the variance tends to infinity. We show on Figure 10 the behaviour of the ratio V ar disc /V ar as a function of µ e and µ a when the mean is kept constant. As a consequence of Chebychev inequality we have (compare with Formulas (24) and (25)):

P(|C -E[C]| ≥ γE[C]) ≤ 1 γ 2 1 1 + r + 2 1 -r 2

Limit when the discretization step tends to 0

The following proposition describes the behaviour of Y (p) when n tends to infinity or, which is the same, when p tends to zero. Note that this also equivalent to letting β e and β a tend to 1.

Proposition 6. Y (p) converges in law to C when p tends to zero.

Proof. The proof is analogous to the one of Proposition 4 in [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF]: one shows that the characteristic function φ p of Y (p) tends pointwise to the characteristic function φ of C. The details are omitted.

In practical terms, this results means that, as long as we are only interested in distributional properties of the concentration, we may consider the discrete model instead of the continuous one provided p is chosen small enough. Note however that the discrete model with arbitrary value of p is interesting in its own right.

Regularity of the discretized model

We assume that µ a ̸ = µ e , since Y (p) is constantly equal to zero when µ a = µ e . The work [START_REF] Albeverio | On fine fractal properties of generalized infinite Bernoulli convolutions[END_REF] provides an analysis of the regularity of generalized infinite Bernoulli convolutions, and specially its fractal properties. The interested reader may also consult [START_REF] Erdós | On a family of symmetric Bernoulli convolutions[END_REF][START_REF] Peres | Sixty years of Bernoulli convolutions[END_REF][START_REF] Peres | Absolute Continuity of Bernoulli Convolution, a Simple Proof[END_REF] for studies on regular infinite Bernoulli convolutions.

The relevant parameter here is β * = β e ∨ β a . The results of [START_REF] Albeverio | On fine fractal properties of generalized infinite Bernoulli convolutions[END_REF] show that the regularity of the discretized model is very different depending on whether β * < 1 2 or β * ≥ 1 2 : • When β * < 1 2 , the distribution of Y (p) is singular with respect to the Lebesgue measure. More precisely, the Hausdorff dimension of its support is equal to lim inf )

. The Hausdorff dimension measures the "size", in a certain sense, of the support of the measure, i.e. where it is concentrated. A Hausdorff dimension smaller than one means that the set of possible values taken by the concentration is extremely sparse and does not form a continuum. In addition, the probability of being in an interval varies in a very non-smooth way with the bounds of the interval. In practical terms, this means that probability distribution of the concentration will be highly irregular : it will vary erratically, taking only very particular values, and the probability that it ranges in some interval may vary wildly when the interval changes. Note that β * < 1 2 is equivalent to p > µ * log(2) or hk * > log(2) with k * = k e ∧ k a , (recall that h is the discretization step ). This reflects the fact that, for fixed k * , the time instants at which the patient is supposed to take his drugs are sufficiently spaced, and that he forgets to do so randomly at some of these instants.

32

Alternatively, for a fixed h, k * must be sufficiently large. Note also that, in view of p < 1, this is possible only if µ * < 1 log(2) . • The case β * > 1 2 is mathematically much more delicate and is not completely settled. Results given in [START_REF] Cooper | Dimension, measure and infinite Bernoulli convolutions[END_REF][START_REF] Pratsiovytyi | A class of random variables of the Jessen-Wintner type[END_REF] entail that the distribution of Y (p) is absolutely continuous for Lebesgue-almost all β * in ( 1 2 , 1) when p = 1 2 . This means that the distribution of the concentration has a probability density: in other words, in contrast with the case β * < 1 2 , we are in the usual situation where the probability that the concentration ranges in some interval varies smoothly when the interval changes.

• In all other cases, the characterization of the regularity of the distribution of Y (p) remains an open problem.

6 Discrete model with random intake instants and random doses

In the discrete case, the random-dose model takes the following form:

• We still have a sequence of i.i.d. Bernoulli random variables (X j ) j≥0 with parameter p, which mark the random instants of drug intake.

• To account for the random dosing, we consider a sequence of i.i.d. random variables (D j ) j≥0 , that will represent the doses, with distribution ν compactly supported on R * + . We will let ν denote the product measure ν ⊗ N .

• We make the assumption E[D 1 ] = 1, which means that the patient takes the normal dose on average.

• We suppose that all the D i are independent of the X j .

• As before, β e = e -keh and β a = e -kah are respectively the elimination rate and the absorption rate for one time step.

We investigate the behaviour of the almost sure limit (with respect to ν) of the steady state concentration Y rd (p) , given by Formula (18) in Section 2.6, i.e.

Y rd (p) = α ∞ ∑ j=0 (β j e -β j a )D j X j .
Independence of (X j ) j∈N and (D j ) j∈N entail that the characteristic function of Y rd (p) is given by ] .

φ rd p (θ) = ∞ ∏ j=0 [ (
As before, we first proceed to investigate the variability in this model, and then the limit when the discretization step tends to zero. We do not characterize the regularity of the concentration, as no results are available in this more complex situation.

Figure 1 :

 1 Figure 1: Eight sample paths of the concentration taking into account variability of clearance and bioavailability, graphed for the first fourteen days. The black solid line is the population mean µpop, the dotted-dashed lines correspond to the confidence bands µpop ± σpop and the dashed lines to the confidence bands µpop ± 2σpop.

Figure 2 :

 2 Figure 2: Eight sample paths of the concentration with Poisson distributed instant of intakes, for the first 14 days. The black solid line is the limit concentration mean E(C), the dotted-dashed lines delineate the confidence bands E(C) ± √ Var(C) and the dashed lines to the confidence bands E(C) ± 2 √ Var(C).

Figure 3 :

 3 Figure 3: Evolution of the variance as a function of µe and µa, when all other parameters are kept constant.

Figure 4 :

 4 Figure 4: Evolution of the variance as a function of µe and µa when mean is kept constant.

Figure 5 :

 5 Figure 5: Average number of doses per hour in the non-compliant case as a function of the number of doses per hour in the fully compliant one yielding same variance for the concentration.

Figure 6 :

 6 Figure 6: Average number of doses per hour in the non-compliant case as a function of the variance of the elimination rate and variance of bioavailability, in a fully compliant population yielding same variance for the concentration.

  µ * = 0.5

Figure 7 :

 7 Figure 7: Histogram of the distribution of C

Figure 8 :

 8 Figure 8: Eight sample paths of the concentration with Poisson distributed instant of intakes and discretely distributed random doses, for the first 14 days. The black solid line is the limit concentration mean E(C rd ), the dotted-dashed lines correspond to the confidence bands E(C rd ) ± √ Var(C rd ) and the dashed lines to the confidence bands E(C rd ) ± 2 √ Var(C rd ).

Figure 9 :

 9 Figure 9: Eight sample paths of the concentration with Poisson distributed instant of intakes and uniformly distributed random doses, for the first 14 days. The black solid line is the limit concentration mean E(C rd ), the dotted-dashed lines correspond to the confidence bands E(C rd ) ± √ Var(C rd ) and the dashed lines to the confidence bands E(C rd ) ± 2 √ Var(C rd ).

Figure 10 :

 10 Figure 10: Ratio V ar disc /V ar as a function of µe and µa when the mean is kept constant.

( 1 -

 1 β e )(1 -β a )(1 + β e β a ) (1 + β e )(1 + β a )(1 -β e β a ). (36)

  1 -p) + p ∫ e iθα(β j e -β j a )u ν(du)

Table 1 :

 1 . PK parameters of imatinib. TV: Typical Value; SD: Standard Deviation; CV: Coefficient of variation.

	Parameter Unit TV	SD	CV (%)
	D	mg 150	0	0%
	k a	h -1 0.61	0	0%
	F		1	0.2462	25%
	Cl	lh -1 14.3 0.3491	36%
	V d	l	347 0.5781	63%

  1 λ = 12h, one obtains for the mean and variance: E pop (E d ) = 0.9576 mg/l, and Var pop (E d ) = 0.1836 mg 2 /l 2 , where we have assumed for simplicity that V d is constant, i.e. Cl = k e V d U ke , where U ke is log-normal with parameters m ke = 0 and σ 2 ke = σ 2

Table 2 :

 2 Numerical values of the parameters of the random doses distributions

( 1 -e -2ket ) -µa 2 ( 1 -e -2kat )] .

This means that the integral of the square of the probability density function is finite.

Acknowledgements

The research of L.J. Fermín was supported by Proyect DIGITEO DIM, ANIFRAC : Uncertaities in processes with fractal characteristics, and by a research grant from the Proyect DIUV N • 2/2011 of the University of Valparaíso.

Variability of the concentration

The mean E rd disc and variance V ar rd disc in the random-dose discretized model are easily seen to be:

Random-dosing thus results in multiplying the variance of the deterministic-dose discretized model by Var(D 1 )

1-p + 1. It is then easy to study the variability of the random-dose discretized model using the analysis developed in Section 5.2 for the discretized model.

Limit when the discretization step tends to 0

The following result can be proved in the same way as Proposition 4 in [START_REF] Véhel | Variability and singularity arising from poor compliance in a pharmacokinetic model I: the multi-IV case[END_REF]: Proposition 7. Y (p) converges in law to C rd when p tends to zero.

Again, this means that, if one is only interested in distributional properties of the concentration, the discrete model may be considered in place of the continuous one provided p is chosen small enough.


Bilal Komati
email: bilal.komati@femto-st.fr

Muhammed R Pac
email: muhammed.pac@mavs.uta.edu

Isura Ranatunga

C Édric Cl

Dan O Popa
email: popa@uta.edu

Philippe Lutz
email: philippe.lutz@femto-st.fr

EXPLICIT FORCE CONTROL V.S. IMPEDANCE CONTROL FOR MICROMANIPULATION

This paper presents a study of different force control schemes for controlling contact during manipulation tasks at the microscale. Explicit force control and impedance control are compared in a contact transition scenario consisting of a compliant microforce sensor mounted on a microrobotic positioner, and a compliant microstructure fabricated using Silicon MEMS. A traditional double mass-spring-damper model of the overall robot is employed to develop the closed-loop force controllers. Specific differences between the two control schemes due to the microscale nature of contact are highlighted in this paper from the experimental results obtained. The limitations and tradeoffs of the two control laws at the microscale due to the presence of backlash are discussed. A simple method to deal with the pull-off force effects specific to the microscale is proposed. Future improvements of the impedance control schemes to include adaptation are discussed in order to handle objects with unknown stiffness.

INTRODUCTION

The development of microsystems (MEMS and MOEMS) has led to increasingly complex microstructures and smaller components. This development has increased the need of microassembly and micromanipulation. In the last few decades, researches have defined many micromanipulation and microassembly techniques [START_REF] Savia | Contact micromanipulation-survey of strategies[END_REF], and have developed many robotic systems for microassembly [START_REF] Dechev | Development of a 6 degree of freedom robotic micromanipulation for use in 3d mems microassembly[END_REF][START_REF] Das | µ 3 : Multiscale, deterministic micro-nano assembly system for construction of on-wafer microrobots[END_REF]. The role of these robotic systems is to perform precise positioning and faster processes.

For automating microassembly or micromanipulation tasks, many of the past work deals with vision-based control [START_REF] Wang | Automated 3-d micrograsping tasks performed by visionbased control[END_REF][START_REF] Anis | Automated selection and placement of single cells using vision-based feedback control[END_REF]. In many cases, vision feedback can only enable position control and ignores the interaction forces like gripping forces and contact forces between the grasped micropart and the substrate. In some other simple cases, force information can be derived from vision information if the contact conditions are simple and the environmental stiffness is known [START_REF] Greminger | Vision-based force measurement[END_REF][START_REF] Liu | Vision-based cellular force measurement using an elastic microfabricated device[END_REF]. However, this becomes difficult if both the gripper and the environment are compliant, and if the stiffness of the microparts is unknown. As a result, new microgrippers with integrated force sensing have been recently developed to enable force-feedback control during manipulation at the microscale, such as FemtoTools microgrippers [8]. Furthermore, at the microscale, interaction forces experience significant nonlinear effects that cannot be simply modeled with a spring, such as adhesive forces. It is notably manifested by pull-off force which can be 84 times the 100µm x 1000µm x 1000µm silicon micropart weight [START_REF] Rabenorosoa | Measurement of pull-off force for planar contact at the microscale[END_REF]. Another important reason to take the microscale forces into consideration is the fragility of the components at small scale. Indeed, microparts and microgripppers can easily break during manipulation and handling.

Issues relevant to robot force control at the macroscale are well known due to extensive research conducted in the 80's and 90's. Several classic control schemes have been developed including explicit force control, admittance control [START_REF] Seraji | Adaptive admittance control: An approach to explicit force control in compliant motion[END_REF], impedance control [START_REF] Hogan | Impedance control -an approach to manipulation. i -theory. ii -implementation. iii -applications[END_REF][START_REF] Singh | An analysis of some fundamental problems in adaptive control of force and impedance behavior: theory and experiments[END_REF], and hybrid force/position control [START_REF] Raibert | Hybrid position/force control of manipulators[END_REF]. Recently, several of these control strategies have been used in microscale [START_REF] Komati | Automated guiding task of a flexible micropart using a two-sensing-finger microgripper[END_REF][START_REF] Lu | A force-feedback control system for microassembly[END_REF] but no comparative study has been yet performed to see whether the physics of the microscale adds anything new to the existing body of knowledge.

This paper investigates differences between force control schemes due to scaling. The work in this paper is limited to studying explicit force and impedance controllers, due to the fact that many other schemes are based on them. For instance, admittance control is similar to the impedance control since the admittance is the inverse of the impedance and the hybrid force/position requires to have many axes and to control some axes in position and the others in force.

The importance of this work is that it provides a good theoretical and experimental base for research using force control for micromanipulation and microassembly applications. The paper is organised as follows. The experimental setup used in this paper and some of the microscale specific effects are explained in Section 2. The system modeling and some force control laws formulation are developed in Section 3. Section 4 presents the experimental results obtained for the force control laws integrated into the experimental setup. Section 5 concludes the article.

FORCE CONTROL AT MICROSCALE

A motivating microscale scenario for this paper is depicted in this section and some of the microscale specific effects are explained in order to take them into consideration for the control design.

Experimental Setup

The hardware shown in Fig. 1 and 2 was used to control the interaction between a MEMS force sensor and a compliant environment. A microstructure made of 100 µm thick Silicon via microfabrication was used as the spring mechanism of the environment. A force sensing probe from FemtoTools AG, Switzerland, used as a manipulator, was attached to a motorized stage and the spring mechanism was attached to a stationary base. The stage was controlled through a National Instruments, Austin, Texas, USA servo motion controller. The force feedback was acquired through an analog input data acquisition board, and the controller was implemented in using NI LabVIEW.

The motorized stage from Newport Corporation, Santa Ana, California, USA, has a minimum incremental motion of 0.1 µm and an uni-directional repeatability of 0.5 µm. The stage also has a backlash of 12 µm.

Scaling Issues

There are many differences between contact models at macroscale and microscale. We can mention the different predominant forces, the different signal to noise ratios for sensors and the added nonlinearities/positioning errors for actuators. Surface force being predominant at the microscale, it is required to evaluate the influence of surface forces during the studied task. These forces appear when a contact between the environment and the force sensor happens. To automatically achieve a contact/non-contact transition at the microscale, pull-off force has to be accounted for. It is the force necessary to break a contact due to stiction and has a predominant role when contact happens between microscale objects.

It was shown in [START_REF] Rabenorosoa | Measurement of pull-off force for planar contact at the microscale[END_REF] that the pull-off force can reach 196µN for a planar 50µm x 50µm silicon contact surface that can typically happen in the present case. During the contact/non contact transition, the breaking of the lateral contact may induce a pull-off force for each side of the contact. In this case, the evolution of the lateral contact force according to the position of the micropart can follow curves in Fig. 3,i.e once a contact (micropart/environment) happens, the pull-off force acts as a sticking effect. As shown in Fig. 3, an opposite force is applied during the withdrawal to break the contact. At point A, the lateral contact force is zero but the contact remains due to adhesive force. The contact is broken at B when enough force is applied to balance the adhesive force.

The signal to noise ratio is another difference between the microscale and macroscale. This is caused by the fact that the signals at the microscale have large bandwidths, small amplitudes and high noise. The noise is due to the type of components used, their sensitivity to the environmental conditions (temperature, humidity, small vibrations) and the microfabrication uncertainties. To deal with these microscale specificities, some filters The bandwidth of microsystems usually contains high frequencies, so using low pass filters will affect the dynamics of the system. In addition, the controller must be able to deal with the noise and must not confuse between contact force and noise.

SYSTEM MODELING AND CONTROL

In this section, the model of the system used in this paper is discussed. Then, explicit force control and impedance control are studied showing the differences between these two strategies.

Model of the system

The contact interaction between a robot and the environment is typically modeled by mass-spring-damper models. In this paper, the proposed models are divided into two parts: an environment part modeled by a mass-spring-damper and a gripper part modeled also by mass-spring-damper. This corresponds to a flexible-gripper, flexible-part interaction scenario typically encountered at the microscale. The proposed model is shown in Fig. 4.

In Fig. 4, m is the mass of the contact surface (gripper + environment), m s , k s and d s are respectively the mass, the stiffness and the damping of the gripper, k e and d e are respectively the stiffness and the damping of the environment, x s is the actual position of the stage which moves the gripper, x d and x are respectively the desired and the actual position of the environment. x, x s and x d are relative positions to the initial positions before contact.

During contact, the dynamics equation is given by (1):

m ẍ = F gripper→m + F environment→m (1)
Where F gripper→m and F environment→m are respectively the forces exerted form the spring and damper of the gripper to the mass and from the spring and damper of the environment to the mass. The dynamics of the system is given by (2):

m ẍ + f m + f e = 0 ( 2
) FIGURE 5. EXPLICIT FORCE CONTROL SCHEME
where f m is the force exerted by the mass on the gripper and f e is the force exerted on the environment. f m given by (3), is the force measured by the force sensor.

f m = -F gripper→m = d s ( ẋ -ẋs ) + k s (x -x s ) (3)
The force exerted on the environment f e is given by [START_REF] Wang | Automated 3-d micrograsping tasks performed by visionbased control[END_REF].

f e = -F environment→m = k e x + d e ẋ (4)
Robots are usually controlled by motor efforts (force or torque). In some applications, and especially at the microscale, the robots are controlled in position or velocity instead. In this case, we change the existing control schemes where the calculated command is the force to a controller in which the calculated command by the controller is an incremental position or velocity of a moving stage. The aim is to control the force exerted to the environment f e in order to guarantee the stability of the grasp in micromanipulation, and also to not break components.

Explicit Force Control

The explicit force control aims to control the force exerted on the environment f e in reference to a desired force f d [START_REF] Volpe | A theoretical and experimental investigation of explicit force control strate-gies for manipulators[END_REF]. Two types of explicit force control exist: force-based and inner position loop based. The most commonly discussed technique in the literature is the force-based explicit force control in which P, PI or PD controllers are used [START_REF] An | Dynamic stability issues in force control of manipulators[END_REF][START_REF] Colgate | An analysis of contact instability in terms of passive physical equivalents[END_REF][START_REF] Eppinger | Understanding bandwidth limitations in robot force control[END_REF]. In this article, the inner position loop based explicit force control is used because the robotic system is controlled by position. The explicit force control scheme used in this paper is given by Fig. 5. The force controller calculates the position command x s using the error ε between the desired force and the measured force exerted on the environment f e and sends it to the controlled stage. If f e is not measured, it should be estimated using the position and another measurement, if possible, f m . The form of the controller changes in function of the application and of the requirements of the system. The most commonly used controllers are a P, PI or PD controllers. These controllers are used for their simplicity and their effectiveness. Other types of controllers could be used as well. An incremental controller like the one used in [START_REF] Komati | Automated guiding task of a flexible micropart using a two-sensing-finger microgripper[END_REF] could also be used (see Fig. 6). The idea of this controller is to compare the desired force to the output force and then the controller sends . The advantage of this type of controller is that it is easy to implement, robust and effective even if the system is nonlinear. The drawbacks of this controller is that the response of the system is slow or it chatters as a function of the gain of the controller (see Section 4). A PD controller could be used if the reference force is constant which is usually the case. In this case, the command could be calculated using [START_REF] Anis | Automated selection and placement of single cells using vision-based feedback control[END_REF].

x s = K p ( f d -f e ) -K d ḟe (5)
Where f e in ( 5) is the measured or the estimated force exerted on the environment.

Impedance Control

The aim of the impedance control is to control the dynamics of contact. The dynamics of contact are set to follow a desired impedance. It has been demonstrated in [START_REF] Volpe | Real and artificial forces in the control of manipulators: Theory and experiments[END_REF] that an impedance controller with force feedback contains an explicit force controller. In the case of impedance control with force feedback, the force is also controlled in reference to a desired force. Using the Laplace transform of (2) and ( 3), ( 6) could be derived. The position x is replaced by X(s) = f m (s) d s s+k s + X s (s) using (3):

ms 2 X s (s) + ms 2 + d s s + k s d s s + k s F m (s) + F e (s) = 0 ( 6
)
where s is the complex argument, X s , F m and F e are the Laplace transforms of x s , f m and f e respectively. The desired impedance proposed by Hogan in [START_REF] Hogan | Impedance control -an approach to manipulation. i -theory. ii -implementation. iii -applications[END_REF], using the Laplace transform, is given by [START_REF] Liu | Vision-based cellular force measurement using an elastic microfabricated device[END_REF].

X d (s) -X(s) = 1 M d s 2 + D d s + K d F e (s) (7)
where X d is the desired position, M d , D d and K d define the desired impedance parameters of the system. Using ( 6) and ( 7), the command sent to the stage x s is given by (8): An external force control loop could be added to the internal impedance controller. The goal of this controller is not only to control the impedance of the contact but also to control the exerted force in reference to a desired force. The final impedance control scheme used in this paper is given by Fig. 7 using the controller equation defined in (8) where the controller has two inputs, the desired position x d , calculated by the force controller, and the force measured by the force sensor f m . The output of the controller x s is the position command sent to the system. The force controller in Fig. 7 could be a P, PD, or PI controller. The proportional force controller is chosen in this case because of its simplicity and to guarantee the stability of the controller which is not the case with the integral controller.

X s (s) = M d s 2 +D d s+K d (M d -m)s 2 +D d s+K d X d (s)- (M d -m)s 2 +(D d -d s )p+(K d -k s ) (d s s+k s )((M d -m)s 2 +D d s+K d ) F m (s) (8)

EXPERIMENTAL RESULTS

In this section, the different types of force control discussed in Section 3 are tested, the dynamical behaviors will be shown and the performances will be compared. These experiments were carried out on the experimental system depicted in Fig. 1 and

Estimation of the force exerted on the environment

Before detailing the experimental results, the force feedback is the force measured by the force sensor f m which is not the same as the force exerted on the environment f e as shown in [START_REF] Dechev | Development of a 6 degree of freedom robotic micromanipulation for use in 3d mems microassembly[END_REF]. Thus, f e should be estimated. Using (2), ( 3) and ( 4) f e could be estimated using two possible equations ( 9) or [START_REF] Seraji | Adaptive admittance control: An approach to explicit force control in compliant motion[END_REF]. In the following two equations, the position x is replaced by X(s) =

In [START_REF] Rabenorosoa | Measurement of pull-off force for planar contact at the microscale[END_REF], f e is estimated using the measurement of the position, x s , the measurement of the stage position given by the encoder, and 9) estimates f e without noise. However, [START_REF] Seraji | Adaptive admittance control: An approach to explicit force control in compliant motion[END_REF] estimates f e using f m and ḟm . Where f m is the force measured by the force sensor which is noisy (for the force sensor used the noise has amplitude of 100 µN). Since ḟm is very noisy, f e will also be very noisy. The amplitude of noise could be reduced, in our case, using a low pass filter to filter the measured force f m . Using the estimation of f e given by ( 9) appears to be a good idea. However, a backlash problem of the stage causes a big problem for estimating f e because the measurement of the stage will not be reliable if many reference signals are applied. Thus, ( 10) is used to estimate f e . Before starting the experiments, the parameters used for control m, d s and d e are calculated theoretically and k e and k s are estimated experimentally. The calculation of the mass was done using the volume of the system. The damping coefficients were calculated using the values of the damping of the Silicon and the sections of the system. k s was estimated by applying a force by the force sensor to a stiff environment and measuring the displacement of the stage.

Experimental results for the explicit force control

As discussed in Section 3, two types of explicit force control are tested and compared: a PD controller and an incremental controller. The response of the system in presence of an incremental control and using the estimation of f e given in ( 9) is shown in Fig. 8.

First, the effect of the step (gain in Fig. 6) of the incremental controller influences the behavior of the system. Big step increments imply a low response time but the amplitude of chattering in the static part is also big. However, small step increments implies a big response time and a smaller amplitude of chattering in the static part. Secondly, in these figures, the effect of backlash appears. The estimated force f e returns to zero while the measured force f m does not. The force feedback being the estimated value of f e using (9), the controller controls the estimated force f e correctly but the estimation of f e is wrong because of the backlash problem.

Using a PD controller could be a solution to have a good response time without chattering in the static part. This type of controller is used because it is simple and gives good performance. The controller used is given by ( 5). The performance is given in Fig. 9. In this case, the estimation of f e given by [START_REF] Seraji | Adaptive admittance control: An approach to explicit force control in compliant motion[END_REF] with a low pass filter is used. Fig. 9 shows that f e and f m follow the same shape and they are almost the same in the static part. However, f e and f m are noisy. The response time of the system is good (≈ 1s). It could be improved by modifying the gains of the PD controller but there is a risk that the system will become unstable in this case.

Experimental results for the impedance control

The impedance control block shown in Fig. 7 and the controller given in (8) are used to control the system. Fig. 10 shows the response of the system using the impedance control scheme defined. The system response is a first order with a response time of ≈ 1s. The behavior is close to the explicit force control. This fact is because the impedance control is used with an external force control feedback. The main difference is in terms of the impedance of the system during contact. The impedance control aims to set the dynamics of the system to a desired impedance defined by the user. Comparing f e in Fig. 9 and Fig. 10, the response of the controlled system is first order in Fig. 10 while in Fig. 9 it exhibits a transient behavior before setting f e to f d .

In impedance control, the desired dynamics of the response of the system is set using three parameters M d , D d and K d . Thus, the controller controls at the same time the impedance of the control and the f e in reference to f d . Some experiments have shown that the use of the proposed impedance control scheme is limited to a flexible environment. However, the stability of the proposed controller is not guaranteed for a very stiff environment as shown in Fig. 11.

A simple comparison between the explicit force control and impedance control results shown in Fig. 9 and Fig. 10 shows that the explicit force controller is simpler than the impedance controller in terms of controller complexity. Both of them are robust. The contact transition for the impedance control is better than the explicit force control. The impedance control is more flexible in design (3 parameters to control) and shows good disturbance rejection.

For instance, the controller is able to set f e to zero when a zero force reference is applied but it is not able to break the contact due to pull-off force. The specific effect due to the microscale is discussed in the following section.

Experimental results showing the pull-off force effect

As shown in Section 2, the pull-off force could appear when a contact happens. The pull-off force acts as a sticking effect. During the experiment, once the contact happens and we aim to separate the contact, setting the desired force to zero will not separate the contact, the contact will persist although the measured force is zero. An easy way to separate the contact is to apply a negative force reference for a short time and then return the desired force to zero. When a negative force reference is applied, the controller aims to set the desired force to the reference and then a command is sent to the stage to come back and then the contact is broken. During this phase, a negative measured force appears due to the pull-off force. This force is not constant during the experiments. It varies as a function of many variables (angle, surface, etc...) [START_REF] Rabenorosoa | Measurement of pull-off force for planar contact at the microscale[END_REF]. Fig. 12, shows the effect of pull-off force on the system. The pull-off force is the negative part of the force f e . It appears while breaking the contact. The contact is not broken directly when the force returns to zero.

Fig. 12 shows also the effect of modifying the estimated values of the parameters. This modification changes the response of the system to a second order with an overshoot in place of a first order response. Other experiments have shown that the impedance controller behavior is acceptable until a parameter change of 40 %. More than 40 %, the reponse of the system will not be acceptable and the system could become unstable.

The problem of signal to noise ratio is avoided, in our case, by adding a low pass filter to the signal because the bandwidth of the system does not contain high frequencies. The results show that the noise level of the output signal f e is small due to filtering.

CONCLUSION

In this paper, a comparison between explicit force control and impedance control for microscale manipulation applications was carried out using experimental results. It was shown that unique effects specific to the microscale affect the effectiveness of these force controllers. Not surprisingly, impedance control showed better performance than explicit force control. Similar to the macroscale, impedance control is better at controlling the dynamics of contact between the manipulator and object. But unlike the macroscale, impedance control is also better at exerting forces on the environment, and more robust at handling flexible components with unknown stiffness. The impedance control scheme used is robust for flexible components and for a small variation of elements. The use of impedance control is limited in presence of a stiff environment. Another limitation would be the parameter uncertainties. Results have shown that the impedance control scheme used in this paper was able to deal with up to 40 % change in the parameters. It was shown that most of the challenges in force control are common for microscale and macroscale but the main differences are the presence of pull-off force at microscale and the signal to noise ratio. A simple method to deal with the pull-off force was proposed. Filtering solves the main problems of the signal to noise ratio, at the expense of bandwidth. To use the impedance control, the parameters of the system should be known m, k e , k s , ... which was the case in this paper. If some parameters are unknown, an adaptive impedance controller could be used to estimate the unknown parameters during the experiments. This will be the focus of future work.

This paper shows that the use of impedance control for microscale is a promising topic in coming years specially with the progress of force control for some microassembly and micromanipulation applications. The impedance control could be used for many microassembly tasks (for picking and releasing a micropart, for guiding a micropart in a rail or for inserting a micropart in a hole). This work was done for a 1 DOF moving stage and it will be extended to more complex systems.

FIGURE 1 .

 1 FIGURE 1. THE DIAGRAM OF THE EXPERIMENTAL SETUP.

FIGURE 2 .

 2 FIGURE 2. SIDE IMAGE OF THE SETUP USED.

FIGURE 3 .

 3 FIGURE 3. THE EVOLUTION OF THE LATERAL CONTACT FORCE IN THE PRESENCE OF PULL-OFF FORCE DURING LEFT SIDE AND RIGHT SIDE CONTACTS IN THE RAIL.

FIGURE 4 .

 4 FIGURE 4. SYSTEM MODEL BASED ON SPRING DAMPER

FIGURE 6 .

 6 FIGURE 6. INCREMENTAL CONTROLLER

FIGURE 7 .

 7 FIGURE 7. IMPEDANCE CONTROL SCHEME

  2 at the University of Texas at Arlington. The control schemes shown in Fig. 5 and 7 consist of motion control, force sensing hardware and a control program developed in LabView as in Fig. 1.

F

  e (s) d e s+k e using (4): F e (s) = d s d e s 2 + (d s k e + d e k s )s + k s k e ms 2 + (d s + d e )s + (k s + k e ) X s (s) (9) F e (s) = d e s + k e ms 2 + d e s + k e F m (s)

FIGURE 8 .

 8 FIGURE 8. SYSTEM RESPONSE FOR AN INCREMENTAL FORCE CONTROLLER (a) WITH A BIG STEP OF THE CON-TROLLER, (b) WITH A SMALL STEP OF THE CONTROLLER

  k e was estimated by applying a force by the force sensor to the flexible environment used and measuring the displacement of the environment. The values used for calculating the control laws are m ≈ 70 mg, d s ≈ d e ≈ 0.3 Ns/m, k e ≈ 20 N/m and k s ≈ 5000 N/m.

FIGURE 9 .FIGURE 10 .

 910 FIGURE 9. SYSTEM RESPONSE FOR A PD FORCE CON-TROLLER

FIGURE 11 .

 11 FIGURE 11. SYSTEM RESPONSE TO IMPEDANCE CONTROL IN PRESENCE OF A STIFF ENVIRONMENT

FIGURE 12 .

 12 FIGURE 12. THE PULL-OFF FORCE APPEARS WHILE BREAK-ING THE CONTACT. EXPERIMENT DONE FOR IMPEDANCE CONTROL SCHEME

ACKNOWLEDGMENT

The authors would like to thank the Net4m project for funding the travel of Bilal Komati from the University of Franche-Comté to the University of Texas at Arlington, the region Franche-Comté in France and the University of Texas at Arlington Research Institute for funding the equipment and hosting the experiments presented in this paper.