Utility of Statistical Model Checking of Stochastic Hybrid Automata for Patient Controlled Analgesia

Opiate concentration in the effect compartment of the brain (OCEC) determines both, the pain control and the side effects. This concentration can be estimated using pharmacodynamics models; however, these models do not predict OCEC when delivery of the drug is random.We are proposing to use stochastic hybrid automata for the verification of individualized model of patient’s drug demands and model of patient’s pharmacokinetics for the estimation of OCEC, and to express results as the probability of falling below the minimum effective analgesic concentration (MEAC) and/or probability of exceeding toxic concentration threshold. Patient controlled analgesia (PCA) model was based on the stochastic hybrid automata, while the verification of the model was done using UPPALL-SMC tool. The suggested approach allowed for quantitative prediction of the OCEC. DOI: http://dx.doi.org/10.5755/j01.eie.23.6.17572

1 Abstract-Opiate concentration in the effect compartment of the brain (OCEC) determines both, the pain control and the side effects.This concentration can be estimated using pharmacodynamics models; however, these models do not predict OCEC when delivery of the drug is random.
We are proposing to use stochastic hybrid automata for the verification of individualized model of patient's drug demands and model of patient's pharmacokinetics for the estimation of OCEC, and to express results as the probability of falling below the minimum effective analgesic concentration (MEAC) and/or probability of exceeding toxic concentration threshold.Patient controlled analgesia (PCA) model was based on the stochastic hybrid automata, while the verification of the model was done using UPPALL-SMC tool.The suggested approach allowed for quantitative prediction of the OCEC.Index Terms-Patient controlled analgesia; hybrid automata; verification; statistical model checking.

I. INTRODUCTION
When patient controlled analgesia (PCA) was originally introduced, the belief was that frequency of analgesic demand uniquely reflects the level of patient's pain [1].Hence, the pain, a very subjective phenomenon, can essentially be formalized "by substitution" of its formalism with the formalism of PCA and by calculating achievable compartmental concentrations.That can be done by describing frequency of demand as a random process.At the same time the prescription can be viewed as an attempt to insure that delivered amount of the drug is sufficient to achieve effective concentration, yet at the same time avoiding side effects [2].The drug delivery with PCA is highly variable and therefore estimation of the probability of drug falling below the minimum effective analgesic concentration (MEAC) and/or exceeding toxic concentration threshold in the effect compartment cannot be done without the model of drug delivery.
There have been efforts made to address safety of the PCA devices utilizing formal verification approach by Sankaranarayanan et al., Arney et al. and Ruksenas et al. [3]- [5].These attempts addressed multiple issues affecting safety of PCA pump, including mechanical plant, operator errors and patient pharmacokinetic differences; however neither of these studies accounted for drug demand variability originating from the individuality of each patient.
Kraujalis addressed the problem of drug delivery and distribution of drug concentration among various compartments using timed automata.However, the main problem with this approach was the necessity to express concentrations and amounts in real numbers: guards and predicates when using time automata formalism are carried out only with integer or rational, therefore real type of data had to be transformed into ranges.Then, using model checking technique for timed automata, the verification could be done using reachability analysis.However, if ranges are broad, an unmanageable "states explosion" occurs.In other words, it can be argued that application of the timed automata for verification of hybrid systems is not an appropriate formalism [6].
We are proposing a novel methodology, which can be used for exactly this purpose, predicting the drug concentration in the effect compartment during variable delivery of the drug.The variable delivery of the drug for individual patient is described using patient's behavior model as statistical ARMA model, which allows to estimate and/predict the stochastic process of drug demand.
Such approach for the first time allowed to combine behavioral model of patient's analgesia demand, pharmacokinetic model (both described as a network of stochastic hybrid automata) with the automatic verification of the model, using UPPALL-SMC, that was developed by the researchers at the Aalborg University and Uppsala University [7], [8].UPPAAL-SMC is a subsystem of UPPAAL which can be effectively used for hybrid systems

Utility of Statistical Model Checking of Stochastic Hybrid Automata for Patient
Controlled Analgesia simulation and verification, using stochastic hybrid automata (SHA) [9], [10].Moreover, the tool UPPAAL-SMC can describe and simulate similar automata models like hybrid automata, price timed automata, or just regular timed automata.The tool also allows checking the properties of those models.Moreover, since the tool UPPAAL-SMC can use real-type variables and does not require transformation of real variables into the integer or rational intervals, it is more versatile than the UPPAAL tool: it allows hybrid models to be accurately modelled and simulated.Moreover, the expansion of UPPAAL tool by adding formalism of stochastic hybrid automata, where clock rates can be arbitrarily changed using various constants and expressions involving other clocks, allows to define ordinary differential equations (ODEs) very effectively [11].The reachability analysis using stochastic hybrid automata is decidable when applying statistical model checking (SMC) technique to the verification of system.The SMC is based on generation of random states trajectories of states graph.Those states trajectories are used to decide whether the system satisfies the property with some degree of confidence.Such application of SMC helps to avoid an exhaustive exploration of the state-space of the system model [8].
Thus UPPAAL-SMC allowed effective verification of our model, opening the possibility of verification of such clinically relevant qualities as the probability of effect concentration falling below the minimum effective analgesic concentration (MEAC) and/or exceeding toxic concentration threshold, all with arbitrary selected confidence intervals [8].

A. Hybrid Automata Model
We used hybrid systems simulation and verification method based on the stochastic hybrid automata (SHA) formalism, which for the clarity purposes we are outlining here.
Hybrid automata (HA) model is a special case of automaton models that is described by a tuple 0 H = ( L,l , X ,E,F ,I ) [12], where  L is a finite set of locations,  0 l is the initial location,  X is a finite set of continuous variables,  E is a finite set of edges where each edge is defined by suite ( l,g ,a, ,l ) a represents an action,  is the binary relation on X R .
 F is the delay function for each l L  ,  I is the invariant mapping that assigns invariant predicates I( l ) to locations, whit the assumptions that: a variable valuation is a mapping v from the continuous variables X to the real's R,v : x R  .So X R represents a set of valuations over X , the delay function F represent evolution of the valuations on X over time.Thus F( d ,v ) is the new valuation on X , after time delay d .Transition between states ( l,v ) of HA is dependent on a continuous value d (delay function) and a discrete value a (action) respectively; where l and v represent valuations for new location and new valuations of variables must exist as an edge e E  , must be available v g  , and ( v,v ) Thus interconnected HA are used for hybrid system modelling and the entire network is called the network of hybrid automata (NHA).NHA is defined as tuple Therefore, when at least one hybrid automaton has stochastic behavior on NHA, then network of HA could be determined as a network of stochastic hybrid automata (NSHA).
The delay density function s  has either uniform or exponential distribution depending on invariant of l .Denote the supremum delay by   D l,v and infimum delay by is an exponential distribution with rate   P l .The output probability function s  is the uniform distribution for every state output probability function s( l,v ) [9].

B. Statistical Model Checking
We used hybrid system verification method, based on statistical model checking (SMC) technique to verify network of hybrid automata models: statistical model checking extends runtime verification capabilities and the result of SMC technique can give answers to whether or not hybrid system satisfies given property or hybrid system satisfies given property with a probability greater than some threshold.Difference between classical model checking and SMC technique is that the result of SMC technique is evaluated with preselected confidence (95 %, 99 %, etc.) [10].
The properties, which are there verified, are described by Weighted Metric Temporal Logic (WMTL).A WMTL formula  over the proposition P and the clock C is generated by the grammar | , where p P  , a b  and c C  .
For these formulas  are obtained from Boolean laws: Formula  means that  should be satisfied starting from the next observation of the run.Formula  is satisfied on the run until 2  is satisfied, and this will continue while value of clock c increases with more than b starting from the beginning of the run, and after it increased for more than a .The release operator R is dual to  , and We used statistical model checking tool UPPAAL-SMC, where model statistical checking technique is realized based on WMTL  which is a subset of WMTL.The logic WMTL  is defined by the grammar where  is satisfied, and this will continue while value of clock c increases with more than d [15].The logic WMTL  is applied for the evaluation of probability or making comparison between them, or to test hypothesis about properties i  .The probability A P ( )  is defined as the probability of the network of stochastic hybrid automata A, that a random run of A satisfies property  . Notation  is satisfied on all runs until clock c is lower than threshold d and the satisfaction of property  is checked [10], [11] and [15].

A. Pharmacokinetic Model
A three compartmental model of drug distribution between the plasma and the brain tissue (effect compartment) was used to describe fentanyl and morphine pharmacokinetics/pharmacodynamics; all of it is shown in Fig. 1.
where 1 x , 2 x , 3 x and e x are the amounts of drug in the central, second, third and effect site compartments, respectively, and 10 k , 12 k , 13 k , 21 k , 31 k and e0 k are the constants defining the elimination as well as intercompartmental transfer rates.The drug delivery function u( t ) is used to represent drug infusion control.
Scheme of patient model analgesia model is presented in Fig. 2.

B. Patient's Behaviour Model
We used two anonymous data logs (both exceeding 24 hour) from the Protocol Library Safety System, Moog, East Aurora, NY.One was with a prescription for morphine and another -for fentanyl.We analysed the length of time period (in minutes) between the two consecutive drug demands [17].Although patient's demands for analgesia may be affected by a multiples factors, they could be reasonably approximated by the autoregressive moving average model of a stochastic process [18].

1) Analysis of Morphine PCA Demands
After comparing 25 different autoregressive moving average models, with parameter values of p and q , ranging from 0 to 4, we chose the model ARMA( p,q ) with lowest value of Bayesian Information Criterion (BIC) is generated [17].The analysis of the length of time period t y (in minutes) between the two consecutive drug requirements was identified, that ARMA(0,0) has the lowest BIC value and it is suggested that the following model is optimal for morphine PCA demands where , identically distributed independent white noise, with mean 0 and variation 2  .
Parameter  is equal to 21,029 and k is 0,407.
Parameter estimates of this model suggest, that lag numerator k does not significantly differ from 0, therefore the time periods t y between two drug requirements during morphine analgesia can be simply modelled as:  , which, according to the value of log likelihood of 222.34, has provided the best fit [17].
Therefore, if r is a basic random number, the noise t  could be generated by the exponential distribution as

2) Analysis of Fentanyl PCA Demands
For the simulation of the time period t y , following model between two consecutive patient demands for fentanyl dose was selected Thus, the value of random variable, having generalized extreme value distribution, can be generated by the following transformation of standard random number r [17]   ln

C. Sensitivity Analysis of ARMA Model
The sensitivity analysis allows to identify model parameters whose inaccurate assessment would significantly reduce the accuracy of the model results.The sensitivity analysis of fentanyl ARMA model ( 7) was performed by varying the ARMA model coefficients and by changing the standard deviation of the noise (9) where ti y is the time period of fentanyl demand estimate by (4) and j ti y -estimate, when varying j -th parameter ( j = 1,2,3 ).Sample size was N = 100 .
Results of sensitivity analysis are showed in Table I, where the first column is variation of the following parameters: coefficients 1 c , 2 c and standard deviation  .
In the columns are presented percentage changes as measured by sensitivity index described by j S (9).The sensitivity analysis demonstrated that small changes in parameter values do not cause extreme changes in the model behaviour.

A. SHA of Pharmacokinetic Model
The results of the implementation of the pharmacodynamics model in the UPPALL-SMC tool, using SHA formalism, are shown in Fig. 3, where ODE system of pharmacokinetic model is described as tool's function pharmODE().Time durations Sup, Mup are estimated as total duration and maximal uninterrupted duration when concentration of drug of the effect compartment is higher than threshold Emax.In order to estimate Mup, an additional variable Tup is used.Variables Sin, Min and Tin represented duration when concentration of the drug in the effect compartment was acceptable, meaning that the drug concentration was higher than MEAC and lower than Emax.The channel inject was used for the connection with pump to initialize virtual drug infusion.We used PBM() from earlier described ARMA model.The duration of drug infusion was described by the variable segment0.Moreover, the PCA model included "lockout interval" which precluded repeat drug injection for I_block minutes.Therefore, if generated duration w was less than I_block minutes, the duration w was changed to I_block minutes.

Developed model of patient controlled analgesia controls infusion of drug and has three states
Drug infusion dose and rate was simulated using: potion/segment parameters (Fig. 3).
On state Sleep the drug administration was not allowed.Action inject! was communicated to the pharmacokinetic model and was set to transmit the signal for the start of drug infusion.

A. Initial Conditions and Assumptions (Simulation Protocol)
Simulation of morphine and fentanyl PCA was performed according to the aggregate scheme [19] as it is shown in the Fig. 1.Drug infusion controller simulated patient's demands, the probabilistic properties of which were defined by the ARMA model that is described in Section II.
Pharmacokinetics of morphine and fentanyl was simulated by the three compartment model.The following parameters were used for the simulation of morphine pharmacokinetics: Central compartment volume = 17.8 l, Time to deliver the bolus dose = 40 sec (infusion rate for Moog, East Aurora, PCA pump), Bolus dose = 1 ml., Lockout interval 10 min.
Morphine micro rate constants were chosen from [20].
For the simulation of fentanyl PCA, the following parameters were used: Central compartment volume = 6.09l, Time to deliver the bolus dose = 40 sec, Bolus dose = 1 ml., Lockout interval 6 min.
Fentanyl micro rate constants were chosen from Shafer publications [16].

B. Model Properties Checking
We checked seven properties of the PCA model for being "true" or "false".Duration of the simulation and checking was chosen to be 24 hours.We denoted u  i.e.: Is the total time more than 90 % checking time (0.9•24 hours), when concentration of the drug in the effect compartment e x is between MEAC and Emax, i.e.:  II.

C. PCA Simulation
Verification time of PCA model depends on the selected simulation/verification period.Some calculations took more than 2 hours when simulation/checking time was 24 hours.
Therefore, we made several steps to reduce computational time required for verification: we simulated ten stochastic runs of PCA model to identify time to reach the steady level of drug concentration in all compartments.
Figure 5 shows one stochastic run of PCA model with fentanyl (using UPPALL-SMC query: simulate 1 [ < = 86400] { SCr.X1, SCr.X2, SCr.X3 } ): green curve is fentanyl concentration in the slow equilibration compartment, yellow curve -drug concentration in the rapidly equilibrating compartment and red curve -in the central compartment.
Ten simulating runs of PCA model demonstrated that drug concentration reaches steady-state after about 2 hours' simulation (checking) time.Therefore, the initial conditions of PCA model could be replaced by the averages of variables at the end of second simulation hour.Modified initial state of PCA model is shown in Table III.After the steady-state of a drug concentration is reached, it fluctuates around constant level.Therefore, the time for simulation and checking of the properties dropped down to 3 hours.The simulation results of ten stochastic runs of PCA model for fentanyl using steady state modification of the initial state (from Table III) are shown in Fig. 6.In Fig. 7 is shown one run of the drug concentration e x for three hours of simulation when modified initial conditions were used.In Fig. 7 Fentanyl concentration in the effect compartment indicates that concentration of drug crosses selected threshold (MEAC or Emax) several times.Besides, the interval, when drug concentration is outside of the interval between MEAC and Emax, is variable.

VI. VERIFICATION RESULTS
We verified all seven PCA model properties using steady state values for the initial state with total checking time of three hours.
Properties     Table V. Morphine PCA satisfies six of seven properties determining with the probability > = 0.78.The third property is not satisfied: it means that the total time period, when concentration of the drug in the effect compartment e x is higher than MEAC, is shorter than 100 %, but longer than 90 % (seventh property shows that 90 % of simulation time, e x is between MEAC and toxic threshold).

VII. DISCUSSION
Earlier research suggests, that different people show different demand patterns when using PCA [1], [21].This behavioural pattern depends on a variety of simultaneously occurring factors: the level of pain, drug concentration in plasma and in the effect site, various side effects (e.g., nausea, sedation, respiratory depression), or even the psychological state of the patient, such as anxiety, cognitive impairment or sleep [2].All these simultaneous factors introduce randomness in the demand pattern [22].Therefore, demand sequence could be thought of as a random process that allows to use auto-regression moving average model ( ARMA( p,q )) which is commonly used in the time series forecasting [18].Depending on the different type of drugs being used, ARMA model exhibits variability of its parameters p and q [17].
Once acquired, ARMA behavioural model can be combined with the hybrid automata pharmacokineticpharmacodynamics model for simulation and estimation of the drug concentration in the effect compartment.This approach allowed predicting of the probability that drug concentration in the effect compartment will fall below the predetermined MEAC (minimal effective analgesic concentration) and/or will exceed Emax (predetermined concentration for the side effects of the drug).
The proposed methodology allows to combine behavioural model of the patient with the formalism of pharmacodynamics model and thus to predict drug concentration variation in the effect-compartment in each individual patient.
Main value of this approach is the ability to formally check the stochastic variations in plasma and effect compartment concentrations in each individual patient against experimentally established populational norms for minimal effective and toxic concentration.If populational norms for the given drug are not available, the proposed approach allows their estimation from the larger sample of PCA logs.
Furthermore, this approach can be used not only for a single drug, but also for the combination of medications, such as morphine/fentanyl as it was used by Friedman [23].When experimentally determining interaction at different ratios is not practical or feasible, the model based checking may be used to find optimal combination.
MEAC and Emax values in this analysis were chosen arbitrarily from the populational data.To determine individual MEAC and Emax (in addition to pump log) we would need more data, namely the information about the level of pain and the degree of respiratory depression.However, even with these limitations, combining pharmacokinetic and patient behavioural model allowed to estimate probabilities and average duration of intervals when drug concentration exceeds populational norms for analgesia and toxic effect thresholds.Described simulation and verification system for the patient controlled analgesia allows probabilistic prediction of the drug concentration that reaches and exceeds predetermined population threshold.This can be used for populational analysis of PCA prescription and to form the basis of a standardized drug library.
In this study we did not perform systematic analysis of ARMA patient demand models and thus did not determine parameter variability between different patients and between different days of PCA therapy.Although it would require larger data set, hybrid simulation method that was used here makes such analysis feasible.
Moreover, the methodology that was used to evaluate patient controlled analgesia performance avoids "explosions of states" that otherwise might have occurred when using time automata verification.In summary, statistical model checking verified the stochastic hybrid automata based analgesia model and allowed statistical prediction of the opiate concentration in the effect compartment.

VIII. CONCLUSIONS
We developed stochastic hybrid automata model and automatic model checking procedure that checks for the concentrations of two drugs in the effect compartment given the personalized pattern of patient demand for these drugs.This can be used for in silico evaluation of pharmacodynamics effects of customary compounded drug mixture, for guiding the transition to oral medication and for drug tapering and rotation.
Modelling results suggest that periods above critical (toxic) concentration threshold of morphine are less frequent, but of a longer duration as compared to fentanyl.
While clinical comparison did not detect difference in the side effects and effectiveness between morphine and fentanyl PCA [24], more frequent dose adjustments were needed for fentanyl PCA.This is consistent with our verification results that demonstrate larger variation in fentanyl concentration.
Special clinical application of this approach allows the investigation infusion of drug mixtures [23] (such as morphine/fentanyl; morphine/alfentanyl; morphine/ketamine for patient controlled analgesia) that can easily cause a variety of the desirable and nefarious side effects simultaneously.
Of note, in some patients MEAC may exceed toxic concentration, leading to analgesia/sedation mismatch, when patient demonstrates toxic effects and inadequate analgesia at the same time.And although individual MEAC and toxic concentration thresholds were not addressed in this study, it is certainly possible to do such estimation.For this purpose, data about pain ratings and side effects, such as respiratory depression, has to be available.
Side effects are random variable which is conditioned on concentration; and although it is true that analgesia and side effects cannot be determined solely by the concentration of the drug, it is also true that without the estimation of the concentration of the drug, no statement can be made about drug effectiveness and the potential for side effects in the particular patient.
Once populational normative for MEAC and toxic concentration are established, model based verification will allow to determine the patients that will need higher monitoring acuity, institution of continuous positive pressure and/or decrease of the PCA dose.As such, the methodology described in this article can be used to standardize PCA prescription using populational data.

2 >=
The notation c d    is used on the statistical model checking algorithms of UPPAAL-SMC and can answer the following three types of questions: Hypothesis testing: Is the probability NSHA A greater or equal to a certain threshold p [0,1] NSHA A ? Probability comparison: Is the probability , on the UPPAAL-SMC tool, followings queries are used for those three types of questions: Pr bound  , where bound defines how to bound the runs [10].
periods between the two drug demands are identically distributed and independent.Analysis of the ARMA model noise t  has shown it to have the exponential distribution with the parameter = / 31.4

Fig. 3 .
Fig. 3. HA of pharmacokinetic model.Hence, Sdown, Mdown and Tdown are the variables of time when concentration of drug in the effect compartment is lower than MEAC.Variables Sin, Min and Tin represented duration when concentration of the drug in the effect compartment was acceptable, meaning that the drug concentration was higher than MEAC and lower than Emax.The channel inject was used for the connection with pump to initialize virtual drug infusion.B.Drug Infusion Model (SHA Model) Initial, Injection and Sleep.Drug injection is made after the demand from the patient.Time duration w represented the duration between two drug demands and was realized by the function PBM().Function PBM() varied depending on the type of drug and individual patient characteristics.Patient controlled analgesia model is shown in Fig. 4.

iT 2 . 3 .
as time intervals when concentration of drug in effect compartment e x is more than Emax.d i T was denoted as time interval when concentration of drug in effect compartment is lower than MEAC.Then we described the checking properties of PCA model: 1.Is the total time less than h = 1  hour, when concentration of drug in effect compartment e x is higher than Emax, i.e.: Is the uninterrupted time less than m = 14.4 minutes (i.e. 1 % of checking time), when concentration of drug in effect compartment e x is higher than Emax, i.e.: Is the total time equal to the checking time (24 hours), when concentration of the drug in the effect compartment e x is higher than MEAC, i.e.:

4 . 1  5 . 6 .
Is the total time less than h = hour, when concentration of the drug in the effect compartment e x is lower than MEAC, i.e.: Is the uninterrupted time less than m = 14.4  minutes, when concentration of the drug in the effect compartment e x is less than MEAC, i.e.: max Is the total time more than 12 hours, when the concentration of drug in the effect compartment e x is between MEAC and Emax?
Fig. 3. Variables Sup and Sdown are described using notations


from the list of chapter 5.2 for fentanyl and morphine are shown in model with the data that was generated randomly and has generalized extreme values distribution with the same, already estimated, parameters.Distribution between residuals from ARMA and data which generated randomly was compared by means of Kolmogorov-Smirnov test -the result failed to reject the null hypothesis.
(2,el residual.Analysis was performed by maximizing the likelihood criterion.Chi square test did not reject our null hypothesis that residuals follow generalized extreme value distribution.Moreover, we compared residuals from fitted ARMA(2,0)

TABLE I .
SENSITIVITY ANALYSIS OF FENTANYL ARMA MODEL.

TABLE III .
INITIAL STATE OF PCA MODEL.
Table IV and in Table V.

TABLE IV .
CHECKING PROPERTIES FOR FENTANYL PCA.

Table
IV indicates that Fentanyl PCA with probability > 0.87 satisfies all seven properties from chapter 5.2.

TABLE V .
CHECKING PROPERTIES FOR MORPHINE PCA.