Optimizing drug infusion schedules towards personalized cancer chronotherapy

Aims Precision medicine requires accurate technologies for drug administration and proper systems pharmacology approaches for patient data analysis. Here, plasma pharmacokinetics (PK) data of the OPTILIV trial in which cancer patients received oxaliplatin, 5-fluorouracil and irinotecan via chronomodulated schedules delivered by an infusion pump into the hepatic artery were mathematically investigated. Methods A pump-to-patient model was designed in order to accurately represent the drug solution dynamics from the pump to the patient blood. It was connected to semi-mechanistic PK models to analyse inter-patient variability in PK parameters. Results Large time delays of up to 1h41 between the actual pump start and the time of drug detection in patient blood was predicted by the model and confirmed by PK data. Sudden delivery spike in the patient artery due to glucose rinse after drug administration accounted for up to 10.7% of the total drug dose. New model-guided delivery profiles were designed to precisely lead to the drug exposure intended by clinicians. Next, the complete mathematical framework achieved a very good fit to individual time-concentration PK profiles and concluded that inter-subject differences in PK parameters was the lowest for irinotecan, intermediate for oxaliplatin and the largest for 5-fluorouracil. Clustering patients according to their PK parameter values revealed two patient subgroups for each drug in which inter-patient variability was largely decreased compared to that in the total population. Conclusions This study provides a complete mathematical framework to optimize drug infusion pumps and inform on inter-patient PK variability, a step towards precise and personalized cancer chronotherapy. Author summary Accuracy and safety of infusion pumps remain a critical issue in the clinics and the development of accurate mathematical models to optimize drug administration though such devices has a key part to play in the advancement of precision medicine. Here, PK data from cancer patient receiving irinotecan, oxaliplatin and 5-fluorouracil into the hepatic artery via an infusion pump was mathematically investigated. A pump-to-patient model was designed and revealed significant inconsistencies between intended drug profiles and actual plasma concentrations. This mathematical model was then used to suggest improved profiles in order to minimise error and optimise delivery. Physiologically-based PK models of the three drugs were then linked to the pump-to-patient model. The whole framework achieved a very good fit to data and allowed quantifying inter-patient variability in PK parameters and linking them to potential clinical biomarkers via patient clustering. The developed methodology improves our understanding of patient-specific drug pharmacokinetics towards personalized drug administration.

Cancer management is challenged by large inter-and intra-patient variabilities in both 2 disease progression and response to treatments. Thus, the quest for accurate and 3 personalized cancer therapies has fostered the development of new technologies enabling 4 multi-type measurements in individual patients and complex drug scheduling. To 5 translate datasets available for an individual patient into personalized therapies and 6 further ensure their precise administration, new mathematical approaches are required. 7 Indeed, systems medicine, that involves the implementation of theoretical approaches in 8 medical research and practice, is critically needed as emphasized in the roadmaps of the 9 Coordinated Action for Systems Medicine (CaSyM) from the European Union 10 (https://www.casym.eu, [1]) and of the Avicenna action (http://avicenna-isct.org/), and 11 in other international consortia [2][3][4][5]. The final aim is a measurable improvement of 12 patient health through systems-based practice which will enable predictive, personalised, 13 participatory and preventive (P4) medicine [6]. 14 Accuracy and safety of infusion pumps are mandatory to ensure that the correct 15 drug dose is delivered to the patient over the intended period. Recurrent incidents 16 related to devices delivering fluids such as nutrients or medications into the body have 17 led the U.S Food and Drug Administration (FDA) to launch in 2010 an initiative to 18 reduce infusion pump risks 19 (https://www.fda.gov/medicaldevices/productsandmedicalprocedures/ 20 generalhospitaldevicesandsupplies/infusionpumps/ucm202501.htm). Many of 21 the reported events are related to deficiencies in the initial design of the device and of 22 the embedded software. Adverse events may also arise from a defect appearing over the 23 device's life cycle due to technical failure or lack of proper maintenance. However, due 24 to the complexity of the equipment, user errors are also common [7]. 25 Optimizing chemotherapeutics index, defined as the ratio between treatment 26 antitumor efficacy and induced toxicities, is complex at multiple levels. First, large 27 inter-patient variabilities are demonstrated in drug pharmacokinetics, tolerability and 28 anti-tumour efficacy [2,[8][9][10]. Next, important intra-patient variabilities arise from the 29 fact that tumour and healthy tissues, rather than being static over time, display 30 time-dependent variations, in particular over the 24h span, which are called circadian 31 rhythms [11]. The circadian timing system controls most physiological functions of the 32 organism resulting in drug Absorption, Distribution, Metabolism and Elimination 33 (ADME) displaying 24h-rhythms with differences of up to several folds between 34 minimum and maximum activities [12,13]. 35 Chronotherapy -that is administering drugs according to the patient's biological 36 rhythms over 24 h-is a growing field in medicine and especially in oncology. Indeed, at 37 least 22 clinical trials involving a total of 1773 patients with different types of 38 metastatic cancers have demonstrated a significant influence of administration timing 39 on the tolerability of 11 commonly-used antitumor drugs [14]. Two randomized phase 40 III clinical trials in 278 metastatic colorectal cancer (mCRC) patients receiving 41 irinotecan, oxaliplatin and 5-fluorouracil showed that cancer chronotherapy achieved an 42 up-to-5-fold decrease in treatment side effects and nearly doubled anti-tumour efficacy 43 compared to conventional administration of the same drug doses [15]. However, a 44 meta-analysis of these two studies combined to another clinical trial involving 564 45 mCRC patients receiving the same drugs (497 men and 345 women in total) concluded 46 that the chronomodulated drug modality significantly increased the efficacy and 47 survival in men while reducing that in women as compared to conventional 48 administration [16]. Such sex-specificity was further validated for irinotecan 49 chronotoxicity in mouse experiments [17] and in a clinical trial involving 199 mCRC 50 patients treated with oxaliplatin (infusion peak 4pm), 5-fluorouracil (infusion peak 4am) 51 and irinotecan given at 6 different circadian times 52 (https://academic.oup.com/annonc/article/28/suppl_5/mdx393.048/4109820).

53
Both studies showed a higher circadian amplitude in females as compared to males and 54 a difference of several hours between the optimal timing of each gender. Furthermore, 55 circadian biomarker monitoring in individual patients recently revealed up to 12 h 56 inter-patient differences regarding the timing of midsleep, the circadian maximum in 57 skin surface temperature or that in physical activity [18]. These investigations have 58 highlighted the need for the individualization of drug combinations and chronoinfusion 59 schemes to further improve treatment outcome, taking into account the patients' sex,  [13]. 70 The Mélodie pump was recently used to infuse those three anticancer drugs directly into 71 the hepatic artery of metastatic cancer patients in the translational European OPTILIV 72 Study [19]. In this study, the plasma pharmacokinetics of oxaliplatin revealed 73 inconsistencies between programmed delivery schedules and observed drug 74 concentration within the patient blood including a delay in the time taken for the drug 75 to be detectable in the blood and unexpected peaks in plasma concentrations during 76 drug infusion. Such inconsistencies between targeted drug exposure patterns and 77 plasma drug levels motivated the design of a mathematical model of fluid dynamics 78 within the pump system presented hereafter. This pump-to-patient model was then 79 connected to semi-physiological PK models to investigate the inter-patient variability in 80 drug PK after hepatic artery administration. Thus, this systems pharmacology study 81 aimed to develop predictive mathematical models allowing for the quantitative and 82 general understanding of i) the pump dynamics, irrespective of the drug delivery device, 83 and ii) patient-specific whole-body PK of irinotecan, oxaliplatin and 5-fluorouracil after 84 drug administration using an infusion pump. Such mathematical techniques would then 85 allow for precise and personalized drug timing.

87
The overall objective of this study was to accurately investigate the inter-patient 88 variability in the plasma PK of the three anticancer drugs administered during the 89 OPTILIV trial. A first strategy consisted in using compartmental PK modelling taking 90 the delivery profiles programmed into the infusion pump as inputs for the plasma 91 compartments. However, such methodology revealed inconsistencies between the best-fit 92 models and the data, including delays of several hours. We then concluded that the fluid 93 dynamics from the pump to the patient had to be quantitatively modelled. Hence, we 94 designed the complete model in two sequential mathematical studies. First, we studied 95 the drug solution dynamics from the pump to the patient blood for which the model 96 was based on partial differential equations. This novel model of the pump delivery 97 system took into account the specificity of the equipment used in order to accurately 98 predict drug delivery in the patients' blood, although it could be easily adapted to any 99 drug delivery devices. Second, we connected this model to compartmental PK models 100 based on ordinary differential equations. This complete framework allowed the 101 investigation of inter-patient variability in drug PK after hepatic artery administration. 102 Pump-to-patient drug solution dynamics The pump-to-patient model is a transport equation representing the dynamics of the 105 drug solution along the administration tube, with respect to time (t) and 106 one-dimensional space (x)(equation 1). x is the distance along the tube from the pump 107 (x = 0) to the patient (x = L). The drug solution was assumed to be incompressible so 108 that the fluid velocity was considered as constant along the whole tube. Thus, the drug 109 concentration in the tube u(x, t) changes with respect to the following equation: with a Dirichlet boundary condition of, where V (t) is the fluid velocity inside the tube, expressed in mm/h. The constant sa = πr 2 is the cross sectional surface area of the tube (in m 2 ), with r being the radius of the tube. The source term S(t) represents the amount of drug delivered according to the infusion profile programmed into the pump and is expressed in mol/h. Initial conditions along the tube are u(x, 0) = 0. The fluid velocity and source terms are controlled by the pump which imposes a fluid delivery rate expressed in ml/h. They are computed by converting the fluid delivery rate into mm/h and mol/h respectively using the tube geometry and the concentration of each drug solution. Hence, model simulations at the end of the tube (x = L) do not depend on the exact geometry of the tube but rather on its total volume. The input function for PK models depending only on quantities at the end of the tube, the original infusion tube which was constituted of two sections of different diameters was simplified in numerical simulations to a tube of radius 1mm and total length 2340mm that had the same total volume as the original set-up. The total tube volume was set to 1.84 mL as in the equipment used in the OPTILIV study. The transport equation with associated initial and boundary conditions can be solved using the classical method of characteristics which gives [20]: where τ (x, t) is the time at which the drug reaching point x at time t initially entered the system i.e. The input function for the PK models corresponds to the rate of drug infusion into the patient (i.e. at x = L) and can be obtained by: Note that, for all drug infusion apart from the glucose flushes, the source term S(t) 111 is proportional to the fluid velocity V (t) as the drug is infused within the tube in the 112 same time as the fluid, so that d(t) is proportional to V (t) once the tube is filled i.e. for 113 times t such that  inside the tube. This remaining drug was flushed out by the glucose rinse subsequent to 130 drug administration which induced a sudden delivery spike in the patient artery 131 (Fig 1c,e,g). The amount of drug in this spike was expressed in percentage of total drug 132 delivered and was estimated to 10.7% for oxaliplatin, 5.36% for 5-fluorouracil and 1.85% 133 for irinotecan.

134
Our systems approach revealed important differences between the intended drug 135 infusion profile and the actual administration into the patient artery. Hence, we 136 developed optimized infusion profiles that strictly achieved the drug administration 137 intended by clinicians. The same equipment was considered to avoid cost of changing. 138 Drug concentrations of the infusion solutions were kept unchanged in order to avoid 139 possible problems of drug stability. In order to administer the drug in the patient's 140 blood following a smooth sinusoidal function, a profile in three parts is required as 141 follows (Fig 2). The first part of the profile is an initial bolus to fill the tube between 142 the pump and the patient with the drug solution. Once the tube is filled, the original 143 sinusoidal profile starts. Then, to solve the problem of the amount of drug left in the 144 tube when the pump stops, the original sinusoidal profile needs to be interrupted when 145 the total drug amount has left the drug bag. Then, a subsequent glucose rinse needs to 146 be infused according to the final segment of the sinusoidal curve in order to deliver the 147 drug remaining in the tube at the correct rate. The pump-to-patient model provided educated predictions of the drug infusion into the 151 patients' blood, which was a prerequisite to study the inter-patient variability in the PK 152 of irinotecan, oxaliplatin and 5-fluorouracil. A compartmental physiological model was 153 designed for each drug and all parameters were fitted for each patient independently.

154
Compartmental models of irinotecan, oxaliplatin and 5-fluorouracil 155 pharmacokinetics 156 PK models represented the drug fate in: the Liver, to accurately represent hepatic 157 delivery, the Blood, the measurement site, and the rest of the body known throughout 158 this paper as Organs. The volume of each compartment was individualised for each 159 patient using Vauthey method for Liver [21], Nadler's formula for Blood [22], and 160 Sendroy method for Organs [23]. Each model assumed that the drug was delivered represented by linear kinetics. Drug clearance included renal elimination for the Blood 164 compartment, intestinal elimination for the Organs compartment and biliary excretion 165 for the Liver compartment. In the absence of quantitative data and to avoid model 166 over-parametrization, circadian rhythms were neglected in the PK models and all 167 parameters were assumed to be constant over the 8-hour time window of PK 168 measurements. Any chemical species bound either to plasma proteins or to DNA was 169 assumed to be unable to move between compartments or to be cleared from the system. 170 For the sake of simplicity, uptake and efflux rate constants were assumed to be equal for 171 Blood-Liver and Blood-Organs transport respectively, for each of the three drugs.

172
Parameter identifiability assessed though sensitivity analysis to cost function variations 173 revealed poor sensitivity of the clearance rate constant in the Organs compartment for 174 the three drugs (cf. Methods). Hence, Organs clearance was assumed to be equal to 175 that of the Liver compartment for irinotecan and oxaliplatin and was neglected for 176 5-fluorouracil [24]. In the model of irinotecan and 5-fluorouracil, poor sensitivity was 177 also obtained for transport parameters between Blood and Organs, so that Organs 178 transport rates were assumed to be proportional to that of the Liver, and the volumes 179 of the compartments were used to scaled parameters. Parameter likelihood profiles 180 analysis revealed that additional constraints were needed to ensure the identifiability of 181 all parameters (see Methods and SI). Hence, information on renal, intestinal and 182 hepatic clearance relative rates was inferred from literature as follows. For irinotecan, 183 CPT11 drug amount though renal clearance and though combined intestinal elimination 184 and biliary clearance were respectively set to 34% and 51% of the total administered 185 dose [25]. As SN38 renal elimination was documented as negligible, the metabolite was 186 considered to only be cleared through Liver or Organs and these cleared amounts were 187 assumed to account for 15% of the total administered dose of irinotecan [25].

188
Oxaliplatin clearance was set such that 55% of the total administered drug amount was 189 cleared via the kidneys [26]. The amount of Pt bound within the Organs or within the 190 Liver was set to 84% and 12% of the total dose, respectively [27]. 5-FU was shown to be 191 mainly cleared through hepatic metabolism, so that the amount of drug cleared though 192 the Liver was assumed to account for approximately 80% of the total dose [24]. 193 The final irinotecan model had six compartments as each of the three Liver, Blood 194 and Organs, had two sub-compartments: the parent drug irinotecan, and its active 195 metabolite SN38 (Fig. 3). Initial irinotecan administered in the liver was assumed to be 196 only in the form of the parent drug. Irinotecan was activated into SN38 via Michaelis

197
Menten kinetics with the parameter estimates K m taken from [28]. SN38 was 198 considered to only be present in its bound form since the bound fraction is reported to 199 be greater than 95% [29]. SN38 clearance terms accounted for SN38 elimination 200 including its deactivation into SN38G though UDP-glycosyltransferases (UGTs) [28]. Semi-physiological model of irinotecan PK. Compartments were minimised to the most important components, Liver to accurately represent drug delivery, Blood which is measurement site and Organs to represent the rest of the body. C i is the rate constant of clearance from compartment i. Irinotecan is bio-activated into its active metabolite SN38. Irinotecan was assumed to be delivered directly into the liver. The oxaliplatin PK model had six compartments corresponding to bound and free 202 platinum (Pt) molecules in the Liver, Blood and other Organs. Oxaliplatin is rapidly 203 metabolised into platinum complex forms [26], which were not distinguished in the 204 current data so as all metabolites of oxaliplatin were assumed to have the same PK 205 properties in the model. Initial oxaliplatin administered in the liver was assumed to be 206 free. Free Pt could bind to proteins and unbinding from proteins was also included in 207 all compartments (Fig 4).

208
The final model for 5-fluorouracil had three compartments. The drug clearance 209 accounted for both drug elimination and drug metabolism in each compartment (Fig 5). 210 Protein binding of 5-fluorouracil was neglected in the model because of the low protein 211 affinity of this drug [30]. Equations for the three models can be seen in SI. Fig 4. Semi-physiological model of oxaliplatin PK. Compartments were minimised to the most important components, Liver to accurately represent drug delivery, Blood which is measurement site and Organs to represent the rest of the body. C i is the rate constant of clearance from compartment i. Each compartment contains a bound and unbound drug fraction and only unbound molecules can migrate between compartments. b and u are respectively the binding and unbinding rate constants of platinum to proteins. Oxaliplatin was assumed to be delivered directly into the liver in its unbound form. and showed a rapid accumulation of both irinotecan and SN38 in the plasma of patients 223 (Fig 6). No obvious impact on irinotecan and SN38 plasma concentrations was observed 224 regarding the time needed to fill the infusion tube or the 30-min glucose delivery spike, 225 as predicted by the pump-to-patient model.  The fit for the oxaliplatin PK model captured all general trends (Fig 7). The model 227 fit for patient 7 did not fully captured the dynamics of total Pt plasma concentration

235
The 5-fluorouracil model showed a very good fit to data, despite a slight systematic 236 under-estimation of the third datapoint in time. It predicted the glucose flush to induce 237 a late spike in plasma drug concentration which could not be seen in the data for all 238 patients, probably because blood sampling frequency was not high enough (Fig 8). This 239 model-predicted spike in 5-fluorouracil concentration changed the tmax value for 240 Patient 5, 6 and 9. The predicted spike AUC was equal to approximately 5% of the 241 total AUC which was in agreement with the pump-to-patient model prediction. This 242 was only calculable for 5-fluorouracil since its elimination was fast enough for its 243 concentration to be close to zero by the time the glucose flush began. The model fit to each individual patient PK data allowed to investigate the 245 inter-patient variability in resulting PK parameters (Fig 9a, b, c). The CV of each PK 246 parameter was calculated among the patient population (see SI). Interestingly, numbers of clusters were determined by minimising the validity index of Fukuyama and 257 Sugeno V F S as described in [31]. Clustering for different numbers of clusters and their 258 respective V F S can be seen in the SI (Supplementary figures 7, 8 and 9). For irinotecan, 259 the minimum value of V F S was achieved for five clusters. One cluster was composed of 260 Patients 1, 3, 5, 7, 8, and 10, the other four patients were in a cluster on their own. The 261 analysis for oxaliplatin concluded to two clusters, a cluster of only one patient, patient 262 7, and the rest of the patients being clustered together. The analysis for 5-fluorouracil 263 revealed four clusters: 5 patients were grouped in the largest cluster (Patients 1, 2, 3, 7, 264 and 10), two patients in the second cluster (Patients 4, 5) and the final two patients 265 were in clusters on their own. Only patients 1, 3 and 10 were consistently clustered 266 together for all three drugs. Once the patient PK parameters had been clustered, the 267 mean of parameter CVs was reassessed for each cluster with 2 or more patients within. 268 Irinotecan mean CV in the largest cluster was 47.18%, which represented a large 269 decrease compared to the mean CV in the entire patient population equal to 73,16%.

270
Oxaliplatin main cluster which was constituted of all patient but patient 7 had a mean 271 CV of 165.58% as compared to 177,87% for the entire population. 5-fluorouracil's 272 largest cluster had a CV of 28.95% and the smaller cluster had a CV of 51.53%, which 273 corresponded to a drastic decrease of inter-patient variability as the population mean 274 CV was equal to 105.69%. All other clusters for each drug had only a single patient and 275 therefore the CV could not be assessed. Inter-patient variability in drug PK parameters. The first line shows parameter variability across the considered patient population for irinotecan (a), oxaliplatin (b) and 5-fluorouracil (c), the colour and symbols represent the clusters each parameter set belongs to. The parameters are named with reference to the schematics of the models, the subscripts refer to the blood (B), organs (O) and liver (L). In the irinotecan parameters, additional subscripts cpt and sn refer to irinotecan and SN38 respectively. The second line shows multidimensional scaling representation of patient clustering based on their PK parameters for irinotecan (d), oxaliplatin (e) and 5-fluorouracil (f), the x refer to the cluster centroids and the points refer to patient PK parameters projected onto 2D plot.

277
Precision and personalized medicine requires accurate technologies for drug 278 administration and proper systems pharmacology approaches for individual patient 279 multidimensional data analysis. Here, plasma PK data of the OPTILIV trial in which 280 patients received irinotecan, oxaliplatin and 5-fluorouracil through a chronomodulated 281 schedule delivered by an infusion pump into the hepatic artery were mathematically 282 analysed. To allow for an accurate analysis of PK patient data, a model of the pump 283 drug delivery was successfully designed and connected to semi-mechanistic PK models. 284 The overall framework achieved a very good fit to individual time-concentration profiles. 285 The validity of the approach was further demonstrated by the improved data fit using 286 the PDE explicit solution connected to PK models compared to PK models directly 287 integrating infusion profiles that were programmed into the pump (see SI). This study 288 gave insights into inter-patient variability and paved the path to treatment optimization. 289 The simulations for the pump-to-patient model showed and quantified a delay irinotecan. Temporal accuracy is key for precision medicine especially in the context of 296 chronotherapy and chronomodulated drug delivery. Thus, the programmation of any 297 drug administration devices need to account for these delays. The pump-to-patient 298 model that we present here allow to adapt any infusion schemes for any drug 299 administration devices in order to properly administer the treatment schedules initially 300 intended by the oncologists.

301
In addition to such "pump-to-body" delay, the increase in free Pt concentration near 302 22:00 shown in the PK data was explained by a spike in oxaliplatin delivery resulting 303 from the glucose rinse flushing out the residual oxaliplatin left within the infusion tube. 304 This phenomenon was well captured and quantified by oxaliplatin PK model which 305 predicted that the quantity of drug delivered in the final spike was equal to 10.7% of the 306 total dose. The model also showed that the t m ax of oxaliplatin plasma concentration 307 was shifted by several hours due to this delivery profile spike. In silico simulations also 308 predicted that the glucose flush would alter the PK of 5-fluorouracil, however the 309 sampling scheme did not cover the time when this would theoretically happen so that physiologically-based PK model of capecitabine, a pro-drug of 5-fluorouracil, was 325 designed for humans [32]. However, the data available in the OPTILIV study would not 326 allow for estimating parameters of such a detailed model. Next, numerous clinical 327 studies have performed compartment analysis of plasma PK data from cancer patients 328 receiving either 5-fluorouracil, oxaliplatin or irinotecan [33]. These models were 329 designed for intravenous injection and could not be readily used for intra-arterial 330 hepatic administration. Thus, the development of new semi-physiological PK models 331 was necessary to include the drug delivery site as a separate compartment, that was 332 different from the Blood compartment for which data was provided. Furthermore, the 333 intention was also to develop more physiologically-relevant models in view of future 334 account of circadian rhythms and vhronotherapy optimization investigations. Indeed, 335 the developed models are called semi-physiological as the compartment volumes 336 together with relative fractions of clearance routes were inferred from literature. These 337 models could then be further extended to physiologically-based models by detailing the 338 "Organ" compartment and be connected to mechanistic PD models to represent 339 organ-specific drug PK-PD towards chrono-administration optimization.

340
Inter-patient differences in maximum plasma drug concentrations and in the time at 341 which it occurred led us to further investigate variability in between subjects. Irinotecan 342 showed the lowest mean variability. Clustering analysis indicated that patients could be 343 classified into five clusters with respect to irinotecan PK parameters. The second largest 344 inter-patient variability was found for 5-fluorouracil. Clustering for 5-fluorouracil 345 showed there was four clusters. Regarding oxaliplatin, there was the largest variability 346 between patients PK model parameters with all parameters showing high variance.

347
Clustering according to oxaliplatin PK parameters split patients into two clusters 348 leading to isolate patient 7. This clustering of the patients led to a reduced inter-patient 349 variability for all drugs, especially for irinotecan and 5-fluorouracil. This decrease in 350 CVs is not unexpected, but the significant level of reduction means this method could 351 be used as a way to stratify patients into treatment groups with less inter-patient 352 variability in PK profiles. The measure of inter-patient variability could be interpreted 353 as indicators of the need for personalisation as high differences between subjects implies 354 high potential benefit of drug administration personalisation. Here, we demonstrated 355 that the PK of all three considered drugs displayed important inter-subject variability. 356 The remaining clinical challenge lays in determining clinical biomarkers for stratifying 357 patients before drug administration, in order to reach the intended plasma PK levels. In 358 order to do so, we performed modelling analyses and identified the critical PK 359 parameters for irinotecan, 5-fluorouracil and oxaliplatin which were the transport 360 parameters between the Blood and either the Liver or the Organs compartments.

362
In conclusion, a mathematical framework was designed to allow for accurate analysis of 363 patient PK data. A model of the dynamics of the drug solution from the pump to the 364 patient's blood was designed, irrespective of the drug delivery device. It was used to 365 represent the chronomodulated drug administration though the Mélodie infusion pump 366 into the patient hepatic artery of irinotecan, oxaliplatin and 5-fluorouracil. The model 367 revealed significant inconsistencies between the drug profiles programmed into the 368 pump which corresponded to the drug exposure intended by clinicians and the actual 369 plasma PK levels. Importantly, it allowed for the design of innovative drug in-fusion 370 profiles to be programmed into the pump to precisely achieve the desired drug delivery 371 into the patient's blood. Next, the pump-to-patient model was connected to 372 semi-physiological models of the PK of irinotecan, oxaliplatin and 5-fluorouracil. The 373 overall framework achieved a very good fit to data and gave insights into inter-patient 374 variability in the PK of each drug. Potential clinical biomarkers for treatment The pharmacokinetic data used in this investigation came from Lévi et al 382 pharmacokinetic investigation [19] and the comparison study companion study of the 383 European OPTILIV trial (ClinicalTrials.gov study ID NCT00852228), which 384 involved nine centres in four countries [34]. The data has been analysed anonymously. 385

397
The superiority of this drug scheduling compared to non-circadian based administration 398 was demonstrated for intravenous administration within several international clinical 399 trials [15]. Between each drug infusion, there was a glucose serum flush which cleared  [19].
Plasma pharmacokinetics (PK) data was gathered after the first dose of irinotecan, 403 oxaliplatin and 5-fluorouracil and measured longitudinally for each individual patient. have a maximum volume of 2 L. The four channels are controlled by four independent 419 mechanisms which control the delivery to the infusion tube (Fig 1). For the OPTILIV 420 study, the infusion tube comprised of two sections, the first was 135mm long with a 421 diameter of 2.5mm, and the second section was 1500mm long with a diameter of 1mm. 422 The two sections had a total volume of 1.84ml. The four pump reservoirs were loaded 423 with irinotecan, oxaliplatin, 5-fluorouracil and 5% glucose solution respectively, with the 424 latter one being used for washes in between drug infusions [35].

Mathematical modelling 426
A pump-to-patient mathematical model was designed as follows, irrespective of the drug 427 delivery device. The drug solutions dynamics from the pump to the patient's blood was 428 modelled using a Partial Differential Equation (PDE) considering time and 1 spatial 429 dimension. This method was chosen as PDEs can take into account both time and 430 space which was key for modelling systems such as pump delivery. The PDE was solved 431 using a backward finite difference method programmed within Python 3.5.2 432 (https://www.python.org/). The drug PK models were based on Ordinary Differential 433 Equations (ODEs) programmed using Python 3.5.2 and solved using the odeint function 434 from the scipy library [36].

435
PK model parameter estimation involved a weighted least square approach, with 436 conditions also placed on the drug clearance routes. The minimization of the least 437 square cost function was performed by the Covariance Matrix adaptation Evolution 438 Strategy (CMAES) within Python which has been shown to be successful at handling 439 complex cost function landscapes [37]. Model goodness of fit was assessed using the sum 440 of square residuals (SSR) and R 2 values. PK model parameter numerical identifiability 441 given the available data was investigated in a two-step process as follows. First, 442 parameter sensitivity regarding the least-square cost function was computed via a global 443 Sobol sensitivity analysis as a necessary condition for identifiability [38]. This method 444 assesses the relative contributions of each parameter to the variance in the cost function 445 obtained when parameter values are varied, and thus allows for the identification of 446 parameters which have no effect on the cost function and are therefore not identifiable 447 from the available dataset. This step allowed a first reduction of the PK models. Next, 448 likelihood profiles of parameters of the reduced models were derived following the 449 procedure outlined in [39]. Additional biological constraints derived from literature were 450 added to ensure numerical identifiability of all parameters. This two-step model design 451 process was undertaken as computing likelihood profiles is associated with a high 452 computational cost.

453
PK models were fit to single-patient plasma PK datasets independently to obtain 454 patient-specific parameter values. Data was available for 10 to 11 patients which was 455 too few to undertake mixed-effect population analysis and to reliably estimate the 456 parameters variance within a patient population [40,41]. Sampling points at 6 hours  Given the relatively small number of patients, the inter-patient variability in parameter values was assessed using a nearly unbiased estimator of coefficient of variation (CV), where µ is the parameter mean, σ the parameter standard deviation and n is the 476 number of patients.

477
Next, fuzzy c-means clustering was used to define patient clusters based on individual PK parameters, for each drug separately. The fuzzy c-means clustering was done using a python library sckit-fuzzy (http://pythonhosted.org/scikit-fuzzy/). The method is based on the determination of cluster centroids and classification of patient parameter vectors into the clusters such that the following quantity is minimised: where nis the number of patients, c is the number of clusters, x i is the parameter vector of patient i, c j is the centroid of cluster j. w ij is the probability of patient i belonging to cluster j and can be expressed as: Note that, for a given patient i, the following holds: wherec is the average of the centroids. The number of clusters were chosen between 2 478 and n-1 inclusively such that the VFS was minimised. Plotting the clustering results 479 was done using a multidimensional scaling (MDS) algorithm which projects 480 multidimensional data onto a 2D plane while keeping distance metric scaled relatively 481 to original data (Python library sklearn.manifold [42]). Correlation coefficients between 482 original Euclidean distance and 2D-Euclidean distance were calculated were high for all 483 models (¿ 0.98) which showed that the MDS projections were accurate [43]. 484 S1 Fig. 1 Patient data best-fit of irinotecan PK model with original delivery profile and not PDE delivery profile. Each subplot represent an individual patient dataset, fit to the model independently. The top figure shows the fit of irinotecan plasma concentration, the bottom figure shows that of SN38, the active metabolite of irinotecan. SN38 data and model simulations include both bound and free SN38. S1 Eq 2. The equations for oxaliplatin PK model. S1 Table 5. Individual Parameter Estimates of oxaliplatin PK model, including mean and CV of each parameter. S1 Table 6. Sum of Square Residuals (SSR) for oxaliplatin PK model, with either the original delivery profile, or that simulated throught the PDE pump-to-patient model. The table also shows improvement in percentages for each patient and average improvement for all patients. S1 Table 7. R 2 values for oxaliplatin model, with either the original delivery profile, or that simulated throught the PDE pump-to-patient model. The table also shows improvement in percentages for each patient and average improvement for all patients. S1 Fig. 2 Patient data best-fit of oxaliplatin PK model with original delivery profile and not PDE delivery profile. Each subplot is an individual patient data, fit to the model independently. The top figure shows plasma ultrafiltrate platinum concentrations, and the bottom figure shows plasma total platinum concentrations. S1 Eq 3. The equations for 5-fluorouracil PK model. S1 Table 8. Individual Parameter Estimates of 5-fluorouracil PK model, including mean and CV of each parameter. S1 Table 9. Sum of Square Residuals (SSR) for 5-fluorouracil PK model, with either the original delivery profile, or that simulated throught the PDE pump-to-patient model. The table also shows improvement in percentages for each patient and average improvement for all patients. S1 Table 10. R 2 values for 5-fluorouracil model, with either the original delivery profile, or that simulated throught the PDE pump-to-patient model. The table also shows improvement in percentages for each patient and average improvement for all patients. S1 Fig. 3 Patient data best-fit of 5-fluorouracil PK model with original delivery profile and not PDE delivery profile. Each subplot is an individual patient data, fit to the model independently. S1 Fig. 4 Parameter Identifiability for irinotecan PK model. S1 Fig. 5 Parameter Identifiability for oxaliplatin PK model. S1 Fig. 6 Parameter Identifiability for 5-fluorouracil PK model.