Benjamin K Schneider

Arnaud Boyer

Joseph Ciccolini

Fabrice Barlesi

Kenneth Wang

Sébastien Benzekry

Jonathan Paul

M Mochel
email: jmochel@iastate.edu

Sebastien Benzekry
email: sebastien.benzekry@inria.fr

DVM Jonathan P Mochel

Ph.D Sebastien Benzekry

Optimal Scheduling of Bevacizumab and Pemetrexed/cisplatin Dosing in Non-Small Cell Lung Cancer

Keywords: Cancers, Optimization, Oncology, Pharmacokinetics-pharmacodynamics, Mathematical modeling Administer, administered, administration, etc Intravenous NSCLC: Non-Small Cell Lung Carcinoma

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Bevacizumab-pemetrexed/cisplatin (BEV-PEM/CIS) combination therapy has been shown to be an effective first line and maintenance therapy for non-small cell lung cancer (NSCLC) in Phase II and Phase III clinical trials [START_REF] Barlesi | Maintenance bevacizumab-pemetrexed after first-line cisplatin-pemetrexed-bevacizumab for advanced nonsquamous nonsmall-cell lung cancer: updated survival analysis of the AVAPERL (MO22089) randomized phase III trial[END_REF][START_REF] Barlesi | Randomized Phase III Trial of Maintenance Bevacizumab With or Without Pemetrexed After First-Line Induction With Bevacizumab, Cisplatin, and Pemetrexed in Advanced Nonsquamous Non-Small-Cell Lung Cancer: AVAPERL (MO22089)[END_REF][START_REF] Galetta | Cisplatin/Pemetrexed Followed by Maintenance Pemetrexed Versus Carboplatin/Paclitaxel/Bevacizumab Followed by Maintenance Bevacizumab in Advanced Nonsquamous Lung Cancer: The GOIM (Gruppo Oncologico Italia Meridionale) ERACLE Phase III Randomized Trial[END_REF]. Pemetrexed inhibits enzymes necessary for pyrimidine and purine synthesis -primarily thymidylate synthase, which is necessary for thymidine synthesis and tumor cell replication [START_REF] Adjei | Pharmacology and Mechanism of Action of Pemetrexed[END_REF]. Cisplatin is an alkylating agent which crosslinks adjacent N7 centers on purine residues, damaging DNA, disrupting repair, and disrupting purine synthesis [START_REF] Poklar | Influence of cisplatin intrastrand crosslinking on the conformation, thermal stability, and energetics of a 20-mer DNA duplex[END_REF][START_REF] Rudd | Persistence of cisplatin-induced DNA interstrand crosslinking in peripheral blood mononuclear cells from elderly and young individuals[END_REF][START_REF] Dasari | Cisplatin in cancer therapy: molecular mechanisms of action[END_REF]. Disrupting DNA substrate supply results in S-phase arrest, DNA repair disruption, and eventually apoptosis [START_REF] Chen | Pemetrexed Induces S-Phase Arrest and Apoptosis via a Deregulated Activation of Akt Signaling Pathway[END_REF][START_REF] Wagner | Cisplatin-Induced DNA Damage Activates Replication Checkpoint Signaling Components that Differentially Affect Tumor Cell Survival[END_REF]. Cisplatin also significantly disrupts calcium and reactive oxygen species (ROS) regulation, inducing cellular lesions which further sensitizes cancer cells to apoptosis [START_REF] Dasari | Cisplatin in cancer therapy: molecular mechanisms of action[END_REF].

In contrast to the effect of PEM/CIS, i.e. DNA damage, bevacizumab is an anti-VEGF (vascular endothelial growth factor) humanized monoclonal antibody. VEGF is an angiogenic potentiator which promotes the growth of endothelial tissue necessary for arteries, veins, and lymphatics.

By limiting neovascular growth, and therefore blood delivery to neoplasms, bevacizumab exhibits limited antiproliferative properties [START_REF] Ferrara | Discovery and development of bevacizumab, an anti-VEGF antibody for treating cancer[END_REF].

More importantly, bevacizumab transiently induces a pruning effect on neovascular beds, which normalizes blood supply to neovascularly dense tissues (i.e. tumors) [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy[END_REF][START_REF] Tong | Vascular Normalization by Vascular Endothelial Growth Factor Receptor 2 Blockade Induces a Pressure Gradient Across the Vasculature and Improves Drug Penetration in Tumors[END_REF][START_REF] Jain | Antiangiogenesis Strategies Revisited: From Starving Tumors to Alleviating Hypoxia[END_REF]. By normalizing blood supply, bevacizumab enhances chemotherapeutic (i.e. PEM/CIS) delivery to neoplasms [START_REF] Kerbel | Antiangiogenic Therapy: A Universal Chemosensitization Strategy for Cancer? Science[END_REF][START_REF] Tolaney | Role of vascular density and normalization in response to neoadjuvant bevacizumab and chemotherapy in breast cancer patients[END_REF]. The effects of BEV-PEM/CIS are generalized i.e. any cell capable of uptaking the drugs are susceptible to their effects, especially rapidly-dividing cells such as myeloid cells [START_REF] Crawford | Myeloid Growth Factors, Version 2.2017, NCCN Clinical Practice Guidelines in Oncology[END_REF].

Accordingly, BEV-PEM/CIS has a narrow therapeutic window and generalized side-effects [START_REF] Galetta | Cisplatin/Pemetrexed Followed by Maintenance Pemetrexed Versus Carboplatin/Paclitaxel/Bevacizumab Followed by Maintenance Bevacizumab in Advanced Nonsquamous Lung Cancer: The GOIM (Gruppo Oncologico Italia Meridionale) ERACLE Phase III Randomized Trial[END_REF].

Previous studies on BEV-PEM/CIS suggest that sequential administration of BEV-PEM/CIS (i.e.

BEV before PEM/CIS) outperforms concomitant scheduling of BEV-PEM/CIS in treating NSCLC [START_REF] Tong | Vascular Normalization by Vascular Endothelial Growth Factor Receptor 2 Blockade Induces a Pressure Gradient Across the Vasculature and Improves Drug Penetration in Tumors[END_REF][START_REF] Imbs | Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma[END_REF][START_REF] Dickson | Bevacizumab-Induced Transient Remodeling of the Vasculature in Neuroblastoma Xenografts Results in Improved Delivery and Efficacy of Systemically Administered Chemotherapy[END_REF][START_REF] Winkler | Kinetics of vascular normalization by VEGFR2 blockade governs brain tumor response to radiation: Role of oxygenation, angiopoietin-1, and matrix metalloproteinases[END_REF]. This makes BEV-PEM/CIS an attractive target for scheduling optimization via modeling and simulation, as a range of practical predictions -such as optimal scheduling in humans -can be made without the considerable time and resource investment required to conduct in vivo experiments. These predictions can be used to guide future studies, greatly accelerating drug development [START_REF] Ciccolini | Improving efficacy of the combination between antiangiogenic and chemotherapy: Time for mathematical modeling support[END_REF].

In our previous work on BEV-PEM/CIS published in Imbs et al. 2017(17), mice with NSCLC tumor xenografts were administered bevacizumab with PEM/CIS combination therapy (CT) in either concomitant, or delayed (i.e. bevacizumab before PEM/CIS) scheduling. The NSCLC tumors had been modified such that tumor growth could be tracked over time via either bioluminescence or fluorescence. Following previous theoretical investigations, the dataset generated from the mice with bioluminescent tumors was used to develop a semi-mechanistic PK/PD model for tumor dynamics in response to BEV-PEM/CIS [START_REF] Benzekry | A new mathematical model for optimizing the combination between antiangiogenic and cytotoxic drugs in oncology[END_REF][START_REF] Lignet | Theoretical investigation of the efficacy of antiangiogenic drugs combined to chemotherapy in xenografted mice[END_REF]. The model was then used to predict the optimal scheduling gap between bevacizumab and PEM/CIS administration.

The aim of this follow-up modeling work was to both refine and expand upon previous results on BEV-PEM/CIS CT using the much larger fluorescence dataset generated in Imbs et al.

2017 [START_REF] Imbs | Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma[END_REF]. We first showed that the semi-mechanistic model previously developed better explained the data than comparable models (i.e. we validated the previously developed structural model). Then, we refined the parameter estimates of the model, and used it to predict the optimal scheduling gap between bevacizumab and PEM/CIS administration. Next, we used stochastic simulations to explore the marginal loss in therapeutic efficacy when BEV-PEM/CIS was administered at a sub-optimal gap, the effect of bevacizumab dose scaling on population optimal gap, as well as the inter-individual variability (IIV) of optimal gap. Lastly, using literature human PK/PD models and parameter estimates, we were able to scale the model to estimate the optimal scheduling of BEV-PEM/CIS in humans.

Methods

Experimental Procedure

Comprehensive details on animals and the experimental procedure are available in Imbs et al. 2017(17). Briefly, on Day 0 of the experiment, tumors (ca 120,000 cells) consisting of H460 human NSCLC transfected with luciferase and the tdTomato gene (H460 Luc+ tdTomato+, Perkin Elmer France) were injected ectopically into the left flank of 90 mice. Animals were pathogen-free, immunocompromised, 6-week-old, female Swiss nude mice (Charles River, France). The mice were randomized into one of five treatment groups. The first study group (Control) received no treatment. The second treatment group (PEM/CIS) was administered both 100 mg/kg of pemetrexed IP and 3 mg/kg of cisplatin IP on Days 14, 28, and 42 of the experiment. The third, fourth, and fifth treatment groups (BEV-PEM/CIS) received the same PEM/CIS treatment as the second experimental group.

In addition, the BEV-PEM/CIS treatment groups were administered 20 mg/kg IP of bevacizumab either concomitantly with the PEM/CIS administrations (Group 3), 3 days prior to each PEM/CIS administration (Group 4), or 8 days prior to each PEM/CIS administration (Group 5) -see Table S1 for administration tabulation.

Tumor growth was monitored on a minimum bi-weekly basis using Ivis Spectrum imager (Perkin Elmer France) and images were acquired and analyzed using the Living Image 6.0. software (Perkin Elmer France).

Mice were supplied with paracetamol supplemented water (e.g. 80 mg/kg/day) to prevent disease-related pain. Animals showing signs of distress, pain, cachexia (i.e. loss of 10% of body mass), or extensive tumor proliferation (i.e. within 2-3 cm) were euthanized. All animals were euthanized on Day 87 of the experiment. All experiments were approved by the local ethical committee at French Ministère de l'Education Nationale, de l'Enseignement Supérieur et de la Recherche, and registered as #2015110616255292.

PK/PD Structural Model Building and Evaluation

The PK models for bevacizumab, pemetrexed, and cisplatin were derived from previously published PK models in mice [START_REF] Johnsson | Pharmacokinetics and tissue distribution of cisplatin in nude mice: platinum levels and cisplatin-DNA adducts[END_REF][START_REF] Lin | Preclinical Pharmacokinetics, Interspecies Scaling, and Tissue Distribution of a Humanized Monoclonal Antibody against Vascular Endothelial Growth Factor[END_REF][START_REF] Woodland | Metabolism and Disposition of the Antifolate LY231514 in Mice and Dogs[END_REF]. The parameters for these models were fixed to the typical values from those studies and assumed no IIV.

The PD model was selected from a series of sequentially fit tumor growth and drug effect models. First -using only the control dataset -the exponential, linear-exponential, and Gompertz growth model were cross-evaluated as models of unperturbed tumor growth [START_REF] Benzekry | Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth[END_REF].

Then, incorporating the full dataset into the fit, the log-kill effect of bevacizumab, log-kill effect of pemetrexed, and log-kill effect of cisplatin were each considered. The interaction effect between bevacizumab and PEM/CIS was included to represent the synergistic effect of bevacizumab.

Following previous work for the effect of cytotoxic drugs, three cellular death compartments were included in the PD portion of the model to represent the delay between cellular damage due to PEM/CIS and cell death [START_REF] Simeoni | Predictive Pharmacokinetic-Pharmacodynamic Modeling of Tumor Growth Kinetics in Xenograft Models after Administration of Anticancer Agents[END_REF].

Competing models were evaluated numerically using Bayesian information criteria (BIC) and the precision of parameter estimates -defined as the relative standard error of the estimate (RSE).

Observed vs. predicted plots, individual fit plots, and Visual Predictive Checks (VPCs) were produced to graphically assist model evaluation (as automated in Monolix 2018R2). VPCs were produced using the default estimation process for VPCs as of Monolix 2018R2 i.e. to create the 90% prediction intervals for the 10 th , 50 th , and 90 th percentiles, 500 simulations are performed using random individual parameters and the design structure of the experiment.

After model selection, the statistical correlations between random-effects were explored via visual inspection. Correlations plots between random effects were produced, defined as 𝜂 ",$,% 1 vs 𝜂 ",$,% 2 i.e. the random effect 𝜂, of individual 𝑖, at time 𝑡, of parameter 1, i.e. 𝜙 1 , vs the random effect 𝜂, of individual 𝑖, at time 𝑡, of parameter 2, i.e. 𝜙 2 . The full posterior distribution of the parameters were used in place of EBEs to avoid visual artifacts due to shrinkage as suggested by Lavielle, et al, 2016(28) and Pelligand L, et al, 2016(29). Statistical correlations between random effects were also numerically assessed using a Pearson correlation test at a P < .05 threshold.

SAEM convergence and final model parameterization were graphically assessed by inspection of search stability, distribution of the individual parameters, distribution of the random effects, individual prediction vs. observation, individual fits, distributions of the weighted residuals, as well as VPCs.

The precision of parameter estimates was numerically assessed using RSE. The normality of random effects distributions, the normality of individual parameter distributions, and the normality of the distribution of residuals were each numerically assessed using a Shapiro-Wilk test (P < .05). The centering of the distribution of residuals (i.e. centered on 0) was numerically assessed using a Van Der Waerden test (P < .05).

Parameter stability was assessed by comparing parameterizations resulting from random-initial starting value selection -as implemented in the Monolix assessment suite. The assessment suite performs 5 SAEM parameterizations in series using random initial parameter values uniformly drawn from the interval from approximately 60% to 160% of final parameter estimates.

The SAEM of the individual parameterizations was then tracked between runs -giving a range of parameter value estimates, RSEs, and log-likelihoods to compare. This assessment was used to ensure that the algorithm did not converge to a log-likelihood local minimum during the process of producing final parameter estimates. The settings described are the default as of Monolix 2018R2.

Simulations

Simulations were performed in R 3.4.4 using Simulx 3.3.0 [START_REF] Lavielle | Simulation of Longitudinal Data[END_REF] to simulate from Monolix run files. First, a function was built which accepted treatment schedule, parameter substitutions, dose, and number of individuals as input and produced a simulated population as an output.

This function was simply a convenience wrapper of Simulx for automation purposes and was verified by reproducing VPCs per treatment group. Simulation set 1 was used to predict the optimal gap between administration of bevacizumab and PEM/CIS. Simulation set 2 produced an estimate of the IIV of the optimal gap. Simulation set 3 examined the anticipated effect of varying the dose of bevacizumab on the optimal gap. Simulation set 4 scaled predictions of BEV-PEM/CIS efficacy to humans. All simulations were of population level response (i.e. simulated without RSE or IIV) except for simulation set 2 (simulated without RSE and with IIV) which used 1000 monte carlo samples. Further details are provided in the supplementary methods.

Quality assessment

All mlxtran and R codes were assessed for quality control by an independent evaluator.

Results

Error Model

Measurement error was best described using a log-normal constant-error model (equation 1). The natural-log of each individual measurement, 𝑙𝑛(𝑦 ". ) with individual i and repetition j, was modeled as a measurement centered on the natural-log mean of 𝑦 ". over j, i.e. 𝑙𝑛(𝑥"), in addition to some residual error, 𝜖 ". , normally distributed, centered on zero, and with standard deviation 𝑎.

𝑦 ". = 𝑥" ⋅ 𝑒 7 89 𝑤ℎ𝑒𝑟𝑒 𝑙𝑛(𝑦 ". ) = 𝑙𝑛(𝑥") + 𝜖 ". | 𝜖 ∼ 𝑁(0, 𝑎) (Equation 1)

The sGOF graphics and numerical analyses supported that a log-constant error model best fit the data. The log-constant error model was not rejected for both the Kolmogorov-Smirnov test for normality (P = .22) nor the Pearson chi-square normality test (P = .0509). In contrast, a constant error model was rejected by these two tests (see Figure S1 for further details).

PK/PD Structural Model Building

No outliers were identified during initial data exploration. Therefore, no collected data were excluded from model building.

The pharmacokinetics of bevacizumab, pemetrexed, and cisplatin were each modeled using one-compartment models with first order IP absorption and first order elimination based on literature descriptions and PK parameter estimates [START_REF] Johnsson | Pharmacokinetics and tissue distribution of cisplatin in nude mice: platinum levels and cisplatin-DNA adducts[END_REF][START_REF] Lin | Preclinical Pharmacokinetics, Interspecies Scaling, and Tissue Distribution of a Humanized Monoclonal Antibody against Vascular Endothelial Growth Factor[END_REF][START_REF] Woodland | Metabolism and Disposition of the Antifolate LY231514 in Mice and Dogs[END_REF]. Random effects (i.e. 𝜂 BC ) were set to 0 as individual PK was not reported in these experiments.

A Gompertz function (Equation 2) was found to best describe the unperturbed tumor growth 𝑉(𝑡), based on its fit performance over competing models, low RSE on parameter estimates, and literature-established descriptive quality.

Eα -β • ln(𝑉 𝑉 L ⁄ )N • 𝑉 (Equation 2)
Due to sparseness in sampling, more complex semi-mechanistic models of growth were not supported by the data. Parameters 𝛼 and 𝛽 represent the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate, respectively. 𝑉 Q is the unit value of relative fluorescence units (RFU) corresponding to one cell, i.e. the proportionality constant between RFU and the number of cancer cells in the fluorescent volume.

𝑉 Q was estimated externally from Monolix by conducting a naive-pooled, linear-regression on the natural-log of the full dataset. The regression gave a rough estimate of 𝑉 R which was then scaled by the approximate number of cells injected at time 0 (ca. 120,000 cells) to derive 𝑉 Q = 5.064 × 10 YZ RFU.

After selecting an appropriate growth model, the log-kill effects of pemetrexed, cisplatin, and bevacizumab were each considered in parallel. The log-kill effect of bevacizumab was estimated as insignificant and removed from the model. The estimation of the log-kill effect of pemetrexed and cisplatin were found to be highly correlated. To reduce model complexity, only their combined concentration, 𝐶(𝑡), and a corresponding log-kill parameter, 𝛾, were considered in the final model (Equation 3). 3) (Equation 4) 𝑦 ". = 𝐹E𝛷 " , 𝑡 ". N • 𝑒 y z 89 , 𝑗 𝜖 {1, … , 𝑛 " }

⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 𝜏 𝑄(𝑡) V Ż e Ż f Ż g Ṅ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 𝜏 klmno • 10 1 + 𝛿 • 𝐴(𝑡 -𝜏) Eα -β • ln(𝑉 𝑉 L ⁄ )N • 𝑉 -γ𝑄𝐶𝑉 𝛾𝑄𝐶𝑉 -𝑘 • 𝑍 e 𝑘 • (𝑍 e -𝑍 f ) 𝑘 • (𝑍 f -𝑍 g ) 𝑉 + 𝑍 e + 𝑍 f + 𝑍 g ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ - - 𝑉(0) = 𝑉 R 𝑍 e (0) = 0 𝑍 f (0) = 0 𝑍 g (0) = 0 - ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ (Equation
𝐴(𝑡)
𝜙 " = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 𝛼 " = 𝛼 BlB • 𝑒 • €,8 𝛽 " = 𝛽 BlB • 𝑒 • •,8 𝛾 " = 𝛾 BlB • 𝑒 • ',8 𝛿 " = 𝛿 BlB • 𝑒 • ƒ,8 𝑙𝑜𝑔𝑖𝑡(𝜏 klmno 8 ) = 𝑙𝑜𝑔𝑖𝑡(𝜏 klmno "…" ) + 𝜂 † ‡…ˆ‰Š ," 𝑉 R " = 𝑉 R BlB • 𝑒 • ‹ OE ,8 𝑘 " = 𝑘 BlB ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , 𝑖 𝜖 {1, … , 𝑁}
Cellular death due to chemotherapeutic treatment was modeled as a three-compartment transition from the growth compartment to death [START_REF] Simeoni | Predictive Pharmacokinetic-Pharmacodynamic Modeling of Tumor Growth Kinetics in Xenograft Models after Administration of Anticancer Agents[END_REF]. The compartments are labeled 𝑍 e , 𝑍 f , and 𝑍 g (numbering respective to their order) and transition between compartments is governed by intercompartmental clearance parameter 𝑘.

𝑘 was not identifiable using the SAEM algorithm. Thus, after a period of manual exploration, 𝑘 was set to the value of 0.3. This choice is consistent with the parameterization made in Imbs et al [START_REF] Imbs | Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma[END_REF]. This choice also limits the total transition time from the tumor mass compartment to cellular death to the order of a day which is consistent with upper limits of cellular death clearance [START_REF] Spencer | Measuring and Modeling Apoptosis in Single Cells[END_REF].

The full tumor size, 𝑁, was the sum of the size of unperturbed cells, 𝑉, as well as the size of damaged cells undergoing cellular death i.e. 𝑍 e + 𝑍 f + 𝑍 g .

No correlations between random effects were statistically significant enough to be included in the final model (Figure S2). Full parameter estimates and model diagram are provided in Table 1 and Figure 1 respectively. Model diagnostics are collected in Figure 2 through Figure 4.

Simulations

Simulating the experimental treatments with a range of administration gaps from 0 to 10 days (step-size = 0.1 day) suggested that the optimal time delay between scheduling bevacizumab and PEM/CIS in mice is 2.0 days (Simulation set 1).

The simulated IIV of the optimal gap was relatively small. Only three values of individual optimal gap were produced. 96.5% of the virtual animals had an individual optimal gap of 2.0 days, 1.0% of the virtual animals had an individual optimal gap of 2.1 days, and 2.5% of virtual animals had an individual gap of 1.9 days (Simulation set 2).

Scaling the dosage of bevacizumab to either 30 mg/kg or 10 mg/kg produced no effect in the estimated optimal gap and produced no effect in the IIV of the optimal gap (Simulation set 3).

Simulations of the typical human response to chemotherapy and bevacizumab were performed using IV administration, two-compartment absorption, and first order elimination models and parameters [START_REF] Latz | Population pharmacokinetic analysis of ten phase II clinical trials of pemetrexed in cancer patients[END_REF][START_REF] Lu | Clinical pharmacokinetics of bevacizumab in patients with solid tumors[END_REF][START_REF] Urien | Pharmacokinetics of platinum after oral or intravenous cisplatin: a phase 1 study in 32 adult patients[END_REF]. Dosage and frequency of administration recommendations for BEV-PEM/CIS were adapted from DailyMed, a product label database maintained by the U.S.

National Library of Medicine (35). Average adult weight and BSA were obtained from Center for Disease Control and literature estimates respectively [START_REF] Mcdowell | Anthropometric Reference Data for Children and Adults: United States[END_REF][START_REF] Sacco | The Average Body Surface Area of Adult Cancer Patients in the UK: A Multicentre Retrospective Study[END_REF]. In constructing that model, we were able to validate and refine previous modeling in BEV-PEM/CIS. Greater precision in parameter estimates was achieved through external estimation of 𝑉 L , external validation of the residual error model, as well as using the larger fluorescence dataset to obtain final parameter estimates. Then, after exploring a range of predictions in mice, we scaled our model to predict optimal scheduling in humans.

The molecular profile of xenografts has a high degree of similarity with the molecular profile of the primary tumors from which the xenograft was derived [START_REF] Wang | Molecular heterogeneity of non-small cell lung carcinoma patient-derived xenografts closely reflect their primary tumors[END_REF] This indicates that there is compositional heterogeneity between xenografts and primary tumors. In addition, human NSCLC H460 cells are an experimentally established paradigm for modeling NSCLC tumors [START_REF] Li | Addition of bevacizumab enhances antitumor activity of erlotinib against non-small cell lung cancer xenografts depending on VEGF expression[END_REF][START_REF] Nukatsuka | Combination therapy using oral S-1 and targeted agents against human tumor xenografts in nude mice[END_REF][START_REF] Xu | Proof-of-Concept Study of Monitoring Cancer Drug Therapy with Cerenkov Luminescence Imaging[END_REF]. Xenografts were grown after subcutaneous implantation in the mice flank, and not directly in the lungs. This is because monitoring tumor growth of lung orthotopic xenografts with non-invasive techniques is limitingly difficult due to both the rapid movements of the chest in mice during imaging, and photon emission attenuation by the high amount of water in lung tissues. Consequently, ectopic xenografts were considered as the more robust experimental model for data collection. Taken together, these considerations motivated our choice of experimental model and are the basis for which we scaled our mathematical model for making predictions in humans.

In the error modeling portion of the experiment, we demonstrated a strategy through which the choice of the error model can be validated externally to the primary dataset by including supplementary data collection in the experimental design. This simplified the error modeling step in the model building process.

The next stage of this study consisted of determining whether the semi-mechanistic model of tumor dynamics in response to BEV-PEM/CIS developed in Imbs et al. 2017(17) best fit the unfit fluorescence data from the same study. During model building, we attempted to balance our model building procedure between model performance (empirical fit) and the underlying biology, an approach often referred to as the middle-out approach [START_REF] Li | A "middle-out" approach to human pharmacokinetic predictions for OATP substrates using physiologically-based pharmacokinetic modeling[END_REF][START_REF] Bon | Mathematical modeling and simulation in animal health. Part III: Using nonlinear mixedeffects to characterize and quantify variability in drug pharmacokinetics[END_REF].

In selecting potential PD models of tumor growth, several semi-mechanistic models were fit to the experimental data. The Gompertz model and linear-exponential model performed comparably. The parameters of the Gompertz model were estimated with greater precision than the parameters of the linear-exponential model (RSE) where the linear-exponential model was fit with a lower BIC than the Gompertz model. Ultimately, the Gompertz model was chosen over the linear exponential model due to the physiological relevance of its construction.

The parameterization of the final model was slightly unstable due to modest overparameterization. To compensate for this, 𝑘 was fixed to a reasonable physiological estimate to improve precision of parameter estimates, and the search for 𝜏 was upper bounded to reduce spurious individual parameter estimates. In addition, the direct anti-proliferative effect of bevacizumab, and individual effects of pemetrexed and cisplatin, respectively, were removed from the model.

The modeling phase of the study resulted in several validations of previous findings. First, we confirmed the validity of the mathematical model previously published in Imbs et al. 2017 [START_REF] Imbs | Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma[END_REF], which we fit to our dataset. In doing so, we reconfirmed the ability of the model to describe BEV-PEM/CIS scheduling. We also reconfirmed the efficacy improvement of BEV-PEM/CIS dosing over PEM/CIS or control. We observed that a 3 day gap in scheduling is superior to both concomitant scheduling and an 8 day gap in scheduling. We were also able to build on previous work by identifying with greater precision the parameters underlying the mathematical model of BEV-PEM/CIS in NSCLC-tumor bearing mice.

In our mouse simulations, the final tumor volume (after 67 days) in the optimal scheduling group with BEV-PEM/CIS (gap = 2.0 days) was 88.5% of the size of final tumor volume in the concomitant scheduling group. This is consistent with our experimental results i.e. that mice administered bevacizumab approximately 2 days before PEM/CIS have a moderately better response (i.e. greater tumor size reduction) to BEV-PEM/CIS than mice who are administered BEV-PEM/CIS concomitantly.

The small magnitude of improvement in 2.0 day gap over concomitant scheduling is likely due to sub-optimal dosage and frequency of administration of BEV-PEM/CIS in the mice. We observed, by exploration, that a more robust preclinical response might be achieved by doubling the frequency of doses in mice and increasing the individual dosages by 50%.

We also found, through simulation, that scaling the dose of bevacizumab had no effect on the optimal gap and that IIV on gap is low.

Predictions made by our model agree with previous findings in BEV-PEM/CIS scheduling. The order of the optimal scheduling delay (2.0 days) is within the 1 to 5 day gap predictions of previous studies [START_REF] Tong | Vascular Normalization by Vascular Endothelial Growth Factor Receptor 2 Blockade Induces a Pressure Gradient Across the Vasculature and Improves Drug Penetration in Tumors[END_REF][START_REF] Dickson | Bevacizumab-Induced Transient Remodeling of the Vasculature in Neuroblastoma Xenografts Results in Improved Delivery and Efficacy of Systemically Administered Chemotherapy[END_REF][START_REF] Winkler | Kinetics of vascular normalization by VEGFR2 blockade governs brain tumor response to radiation: Role of oxygenation, angiopoietin-1, and matrix metalloproteinases[END_REF]. Studies in tumor perfusion and bevacizumab showed Day 1 and Day 4 decrease in tumor perfusion which is consistent with the marginal predictions in our study i.e. optimal perfusion should be on the order of 2 days with comparable marginal losses on either side of that minimum [START_REF] Van Der Veldt | Rapid Decrease in Delivery of Chemotherapy to Tumors after Anti-VEGF Therapy: Implications for Scheduling of Anti-Angiogenic Drugs[END_REF].

After exploring various predictions in mice made by the model, the PK portion of the model was re-parameterized and the parameters of the PD portion of the model were scaled to simulate the relationship between varied administration schedules (i.e. gap) and efficacy in humans. Using this adapted parameterization, we estimated both optimal schedule of administration of BEV-PEM/CIS in humans and the marginal effects of a sub-optimal administration schedule of BEV-PEM/CIS in humans (i.e. bevacizumab administered too many or too few days before PEM/CIS).

In our human simulations, we predicted a robust improvement in response to sequential BEV-PEM/CIS relative to concomitant scheduling. The final tumor volume (after 67 days) in the optimal scheduling group (gap = 1.2 days) was 30% of the size of the final tumor volume in the concomitant scheduling group. If these predictions are accurate, scheduling optimization could result in significant improvement in BEV-PEM/CIS CT efficacy with no increase in toxicity.

The predicted scale of the increase in efficacy after scheduling optimization was much greater in the human simulations than was either predicted by the mouse simulations, or was measured empirically in mice. This is at least partially due to the greater study that has been undertaken in optimizing dosages for human treatment of NSCLC. However, these promising estimates should not be substituted for clinical testing.

When exploring marginal efficacy loss in sub-optimal administration schedules, we consistently found that the marginal cost of scheduling bevacizumab and PEM/CIS too close together in time was greater than the marginal cost of scheduling bevacizumab and PEM/CIS with too great of a gap in administration -in both mice and humans. This indicates that any potential clinical studies in antiangiogenics and cytotoxics should weight scheduling recommendations toward scheduling at slightly too large of gap.

Finally, the tumor microenvironment is known to be complex and varied. Tumor tissues contain necrotic pockets, heterogenous and dynamic microvasculature, and various sub-mutations which result in differential local growth rate and drug sensitivity. Considering this biological heterogeneity would greatly improve future model predictions and scalability between species.

In summary, our analysis confirms previous findings in BEV-PEM/CIS scheduling while improving precision of parameter estimates, improving prediction quality and detail, and scaling the model to predict the optimal scheduling of BEV-PEM/CIS in humans.

Antiangiogenics will continue to be useful agents in oncology. There are currently several other antiangiogenics regularly used in combination with cytotoxics which could potentially benefit from sequential administration (i.e. antiangiogenic then cytotoxic) [START_REF] Gougis | Clinical pharmacology of anti-angiogenic drugs in oncology[END_REF]. Of note, bevacizumab is currently only approved for concomitant administration with chemotherapy in all of its indications e.g. lung cancer, breast cancer, gastric cancer, etc. This contrasts with the optimized sequential scheduling that model simulations suggest.

There is a recent trend to develop model-informed drug development to optimize anticancer therapy. Our work highlights how mathematical modeling could help to refine clinical treatment modalities. The semi-mechanistic nature of this model allows it to be modularly reconfigured to extend predictions to other antiangiogenics as well as novel therapeutic paradigms such as the immune checkpoint inhibitor, monoclonal antibody pembrolizumab [START_REF] Pai-Scherf | FDA Approval Summary: Pembrolizumab for Treatment of Metastatic Non-Small Cell Lung Cancer: First-Line Therapy and Beyond[END_REF]. This work continues to lay the foundation for building systems pharmacology models of the effect of antiangiogenic and antiproliferative combination therapy in advanced NSCLC. Tortuous vasculature is a phenotype exhibited by many solid tumors, and predicting optimal antiangiogenic scheduling could greatly increase the efficacy of future oncology therapeutics and combination therapies [START_REF] Reck | Pembrolizumab versus Chemotherapy for PD-L1-Positive Non-Small-Cell Lung Cancer[END_REF].

  represents the plasma concentration of bevacizumab. 𝑄(𝑡) represents the synergistic effect of improved vascular quality. In brief, the increase in neoplasm vascular quality due to bevacizumab typically occurs within a period of a few days after administration. To represent this delay in effect, time (𝑡) was delayed by 𝜏. Parameter 𝛿 represents the proportional increase in PEM/CIS efficacy due to vascular quality improvement under bevacizumab therapy.The estimation of 𝜏 was bounded between 0 and 10 using the link function 𝜏 = 𝜏 klmno ⋅ 10, where 𝑙𝑜𝑔𝑖𝑡(𝜏 klmno ) ∼ 𝑁(0,1). All other parameters were best estimated as log-normally distributed. The full statistical representation of individual parameters, 𝜙 " , estimated via SAEM is shown in Equation4, where the full structural model is denoted by 𝐹.

Figure 2

 2 Figure 2 Standard Goodness-of-Fit Diagnostic Plots On the left is individual predictions vs. observations and on the right are the individualized weighted residuals (IWRES) vs time. During model fitting, observations were natural-log-transformed to stabilize predictions. Therefore, residuals, predictions, and observations are natural-log-transformed in these figures. The predictions are approximately normally distributed. On the left, the one-to-one prediction line is the center solid black line, the spline (average agreement between individual prediction and observation) is solid orange, the dashed black lines are the borders of the 90% prediction interval. On the right, the zero residual error line is the center dashed black line, the spline is solid orange. The dashed black lines are the borders of the 90% prediction interval.

Figure 3 Figure 4

 34 Figure 3 Sample of Individual Fits The blue dots represent individual observations while the solid violet line represents individual fits.

Figure 5

 5 Figure 5 Human Pharmacodynamic and Pharmacodynamic Simulations SummaryTo produce these figures, bevacizumab (BEV) was administered anywhere from 0 to 10 days (in steps of 0.1) before pemetrexed/cisplatin (PEM/CIS) was administered. Tumor growth was simulated from 0 to 67 days with no IIV and no RSE. In the top figure, AUC of tumor growth vs gap (0 to 10 days) is depicted. In the middle figure, tumor dynamics over time, with gap indicated by color, are depicted. In the bottom panel, the PK of BEV-PEM/CIS is depicted with gap indicated by color. The top figure indicates that the optimal scheduling gap is 1.2 days and the middle figure depicts the difference in tumor volume between administration gaps. The patient with optimal scheduling had a final tumor volume approx. 30% the size of the concomitant scheduling.

Table 2

 2 Except for the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate (i.e. 𝛼 and 𝛽), the PD model and parameterization were reused exactly as they were determined in the mouse portion of the model. 𝛼 and 𝛽 estimates were obtained from

	Bilous et al. 2018(38), where clinical NSCLC doubling times reported in Friberg and Mattson,
	1998(39) were used to estimate population 𝛼 and 𝛽 for NSCLC in humans. The value of VC
	came from the classical assumption that a 1 mm 3 volume of tumor cells is approximately 10 6
	cells(40). V0 was arbitrarily set to 3 cm 3 .
	The full PK/PD model was then used to simulate the typical cancer growth under various
	administration schedules with a starting tumor volume of 3 cm 3 . Parameter estimates are
	reported in

and simulation summaries are depicted in Figure

5

. The estimated optimal gap between bevacizumab and PEM/CIS administration in humans was 1.2 days.

Discussion

By normalizing tumor vasculature, bevacizumab improves delivery of PEM/CIS to tumors which increases PEM/CIS anti-tumor efficacy. Pemetrexed and cisplatin each have a narrow therapeutic window, low clearance, and high toxicity. It is therefore critical that BEV-PEM/CIS doses are administered as efficiently as possible. This makes BEV-PEM/CIS a natural fit for modeling and simulation studies, as the drug scheduling can be optimized without the need for multiple time and resource intensive in vivo studies. In this analysis we conducted an in silico study of the optimal administration of BEV-PEM/CIS in a xenograft, and human, model of NSCLC by constructing a mathematical model of tumor dynamics in response to BEV-PEM/CIS.

Table 1

 1 Pharmacodynamic Model Parameters for Tumor Proliferation in NSCLC-Xenografted Mice Model parameter estimates (fixed and random-effects) as well as standard errors as determined by the Stochastic Approximation Expectation-Maximization algorithm as implemented in Monolix 2018R2.

	Model Parameter	Estimate			SE	RSE (%)	IIV	IIV (CV%)
	𝜶	0.77 day -1			0.0081	1.06	0.037	4.76
	𝜷	0.04 day -1			0.0005	1.16	0.015	33.55
	𝜸	35.74 (mg•day)	-1	6.36	17.79	0.46	1.28
	𝜹	3.73 (mg•day)	-1	0.92	24.70	0.35	9.51
	𝝉 bound	.119 days			0.0017	1.50	0.017	14.50
	𝑽0	3.4 RFU*			0.27	22.55	1.73	143.03
	𝒌	0.30 days			-	-	-	-
	Residual Error Variance	0.43			0.01	2.55	-	-

Table 2

 2 Model Parameterization for Simulations of NSCLC treated with BEV-PEM/CIS in Humans Except for the proliferation rate of the tumor cells and rate of exponential decrease of the tumor relative growth rate (i.e. 𝛼, 𝛽), the PD model and parameterization were reused exactly as they were determined in the mouse portion of the model. 𝛼 and 𝛽 estimates were obtained from Bilous et al. 2018[START_REF] Bilous | Computational modeling reveals dynamics of brain metastasis in non-small cell lung cancer and provides a tool for personalized therapy[END_REF], where clinical NSCLC doubling times reported inFriberg and Mattson, 1998(39) were used to estimate population 𝛼 and 𝛽 for NSCLC in humans. The value of 𝑉 𝑐 came from the classical assumption that a 1 mm 3 volume of tumor cells is approximately 106 cells[START_REF] Spratt | Rates of growth of human solid neoplasms: Part I[END_REF]. 𝑉 0 was arbitrarily set to 3 cm 3 . 𝑉, 𝑘 𝑚,𝑛 , and 𝑄 represent volume, compartmental clearance from compartment m to compartment n, and intercompartmental clearance respectively.

				𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒		𝛼
				𝑉			𝛽
				𝑄			𝛾
				𝑉			𝛿
							𝜏
				𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒		𝑉
				𝑉			𝑘
				𝑄			𝑉 𝑐
				𝑉
	Pemetrexed Pharmacokinetics		Bevacizumab Pharmacokinetics	Pharmacodynamics
	𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒	91.6 mL/min	𝛼 𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒 𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒	0.22 L/day	𝛼	0.0284 day -1
	𝑉1	12.9 L	𝑉	𝑉1	𝛽	2.8 L	𝛽	1.03E-3 day -1
	𝑄	14.4 mL/min	𝑄	𝑘12	𝛾	0.22 day -1	𝛾	35.7 (mg•day)	-1
	𝑉2	3.38 L	𝑉	𝑘 21	𝛿	0.22 day -1	𝛿	3.73 (mg•day)	-1
	Cisplatin Pharmacokinetics			𝜏		𝜏	1.19 days
	𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒	37.0 L/h	𝐶𝑙𝑒𝑎𝑟𝑎𝑛𝑐𝑒	𝑉		𝑉 0	3 cm	3
	𝑉1	22.7 L	𝑉		𝑘		𝑘	0.3 days
	𝑄	10.7 L/h	𝑄		𝑉 𝑐		𝑉 𝑐	10E-3 mm
	𝑉2	19.5 L

Table S1

 S1 Administration (adm) schedule for the five experimental groups. The cytotoxics consisted of an administration of both 100 mg/kg of pemetrexed IP and 3 mg/kg of cisplatin IP. Bevacizumab was administered in doses of 20 mg/kg IP.

	Group	Experimental Day of	Experimental Day of	Group Size
		bevacizumab adm	cytotoxics adm
	Control	-	-	15
	Cytotoxic control	-	14, 28, 42	15
	No gap in adm	14, 28, 42	14, 28, 42	16
	Three day gap	11, 25, 39	14, 28, 42	15
	Eight day gap	6, 20, 34	14, 28, 42

*RFU: Relative Fluorescence Unit

•cell -1 𝑉

Acknowledgements

We would like to thank Yeon-Jung Seo of Iowa State University College of Veterinary Medicine, Ames, IA, U.S.A, for performing the quality control assessment on code written to produce this study.

We thank Dr Valerie Le Morvan and Professor Jacques Robert from the Institut Bergonié, France, for their kind assistance.

The work reported in this manuscripted was not funded or sponsered.

Table S2 Administration (adm) schedule for the simulated human individual. Drugs were administered on three separate occasions, the 17 th , 31 st , and 45 th experimental day. Bevacizumab administration (adm) was varied between simulation run i.e. it was administered anywhere between 0 and 10 days before Pemetrexed-Cisplatin.

Bevacizumab

Pemetrexed Cisplatin

Adm Day 1

Figure S3

Below are various SAEM convergence results produced by Monolix's assessment suite. The best performing (defined as lowest -2 loglikelihood) parameterization (in yellow) closely matches our final parameter estimates. The assessment indicates that there is some degree of instability in parameter estimates. But, variance in parameter estimates is relatively low and RSE on final parameter estimates is low, indicating that we've achieved a reasonable balance between the empirical fit and underlying biology.

Especially important to our analysis is the relatively low variability in Tau parameter estimates.

Figure S4

Below is the 3D surface of simulated tumor growth response (dimension 1) with respect to gap (dimension 2) and time (dimension 3) in mice. This 3D surface allowed us to more intuitively understand the effect of varied administration schedules. It is available for manipulation at https://plot.ly/~Benjamin-PKPD/15/#/.

Figure S5

Below is the 3D surface of simulated tumor growth response (dimension 1) with respect to gap (dimension 2) and time (dimension 3) in humans. This 3D surface allowed us to more intuitively understand the effect of varied administration schedules. It is available for manipulation at https://plot.ly/~Benjamin-PKPD/17/#/.

Supplementary Methods: Simulation Details

In simulation set 1, population response was simulated with no IIV and without RSE. The simulated treatment schedule consisted of 100 mg/kg of pemetrexed IP and 3 mg/kg of cisplatin IP on Day 14, 28, and 42. 20 mg/kg IP of bevacizumab was virtually administered anywhere from 0 to 10 days (in steps of .1) before Day 14, 28, and 42. See Table S2 for details.

The simulation was stopped after 67 virtual days. These simulations allowed the prediction of the optimal gap between administering bevacizumab and PEM-CIS. Efficacy of the drugs was assessed using AUC of the tumor growth over time as calculated by the spline method implemented in the R package DescTools 0.99.25 1 . The optimal gap was defined as the administration scheduling gap which resulted in the lowest AUC of tumor growth vs. time over the simulated time period (t = 0 to t = 67).

In simulation set 2 the inter-individual variability of optimal gap was meticulously calculated using a virtual population of 1000 simulated mice. First, a set of parameters representing 1000 virtual mice was generated without RSE from the parameter distributions determined in the structural modeling phase. Then, the procedure from simulation set 1 was replicated per individual to estimate each individual optimal gap. Finally, the distribution of individual optimal gaps was assessed using the mode, standard deviation, and range as well as graphical methods such as a violin plot of count density.

To produce simulation set 3, the procedure which generated simulation set 1 was repeated with scaled doses of bevacizumab. First, 10 mg/kg IP of bevacizumab was substituted for 20 mg/kg IP of bevacizumab. Second, 30 mg/kg IP of bevacizumab substituted for 20 mg/kg IP of bevacizumab. This was done to assess the effect of bevacizumab dosage on the estimated optimal gap.

In simulation set 4 an analogous treatment regimen to simulation set 1 was simulated in a virtual NSCLC tumor-bearing human population using adapted human dosing recommendations. The random effects were set to zero and population response was simulated Tumor dynamics in response to PEM-CIS-BEV were simulated using a scheduling gap of 0 days to 10 days in steps of .1. A graphical representation of the AUC vs. administration gap was generated to determine the optimal gap between bevacizumab and PEM-CIS dosing.


ProSeg: A comparative corpus of spoken L2 French

Elisabeth DELAIS-ROUSSARIE 1 , Tanja KUPISCH 2 , Paolo MAIRANO 3 , Fabian SANTIAGO 4 & Frida SPLENDIDO 5

1 CNRS & Université de Nantes, 2 Universität Konstanz/UiT The Arctic University of Norway, 3 Université de Lille, 4 Université Paris 8, 5 Lund University . CONTACT: elisabeth.delais-roussarie@univ.nantes.fr

Main features of ProSeg

General features

The

Introduction

Advantages of speech corpora in L2 acquisition studies [2] • They represent of the oral dimension in an L2 in an ecological way (i.e. with spontaneous speech);

• They put in perspective several explanatory factors affecting the L2 acquisition process and competence, such as L1 transfer, speech style (reading vs. spontaneous speech), age of acquisition. [5, 6] • They can be used for cross-comparisons between language pairs, especially when using the same recording protocol [6] • Given these advantages, many current studies are based on speech corpora. But different procedures may be used for data collection, each having its advantages [3] • Experimental data targeting specific linguistic phenomena;

• Naturalistic speech from an oral corpus. ProSeg (Prosody Segmentals) combines the two procedures.

Studies in progess/preparation

(Non-native) gemination in L2 French Foreign Accent Rating

Prosodic phrasing and rhythm

• Prosodic analysis of the segmentation, durational and tonal patterns in accentual phrases (APs) • Data taken from the long text and the monologues.

Comparison of L1 German, L1 Swedish, L1 Italian and L1 French speakers. • Hypothesis: L1 Swedish learners have advantages in producing e APs in L2 French due to the existence of a tonal accent in Swedish.

P ic t u r e d e s c r ip , o n t a s k

  ProSeg Corpus has been designed to allow for : • prosodic and segmental descriptions of L2 French • comparison of potential effects of L1 transfer • contrastive analyses of the oral productions in L1 and L2 with a comparable set of data • analyses of task effects Procedure 1. Background questionnaire 2. Placement test (DIALANG, yes-no vocabulary task, [1]) 3. Reading tasks in L2 French (several texts including a dialogue) 4. Picture description task 5. Free speech (summary: book/ movie) 6. Reading task in the L1 Participants • 65 participants (so far) • Recorded both in L2 French and in their L1 • All learners are university students REFERENCES [1] Alderson, J. C. 2005. Diagnosing Foreign Language Proficiency: The Interface between Learning and Assessment. London: Continuum. [2] Hawkins J.A. & P. Buttery. (2009). Using learner language from corpora to profile levels of proficiency. Studies in Language Testing. Cambridge University Press. [3] Delais-Roussarie E. & H. Yoo (2014). Corpus and research in phonetics and phonology (I): Methodological and formal considerations. In Durand J., Gut U. & G. Kristoffersen (eds), Handbook of Corpus Phonology. Oxford University Press. [4] Mairano, P., Santiago, F. & Delais-Roussarie, E. (2018). Gémination non-native en français d'apprenants italophones. Actes des Journées d'étude sur la parole (JEP 2018), Aix en Provence, Juin 2018 [5] Santiago, F. & E. Delais-Roussarie. (2015a). The acquisition of question intonation by Mexican Spanish Learners of French. In Delais-Roussarie, E., M. Avanzi & S. Herment (eds.), Prosody and languages in contact: L2 acquisition, attrition, languages in multilingual situations. Springer Verlag. [6] Santiago, F. & E. Delais-Roussarie. (2015b). What motivates extra-rising patterns in L2 French: Acquisition Factors or L1 Transfer?, Proceedings of the 18th International Congress of Phonetic Sciences (ICPhS).


Assyr Abdulle
email: assyr.abdulle@epfl.ch

Grigorios A Pavliotis
email: g.pavliotis@imperial.ac.uk

Gilles Vilmart
email: gilles.vilmart@unige.ch

Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations

Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations

Introduction

Sampling from probability measures in high dimensional spaces is an important problem that arises in several applications, including computational statistical physics [START_REF] Lelièvre | Free energy computations[END_REF], Bayesian inference [START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF], and machine learning [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. Typically one is interested in calculating integrals of the form πpφq :" E π φ :"

ż R d φpxq πpdxq, (1)
where πpdxq " πpxq dx 4 is a probability measure in R d , known up to the normalization constant and φ P L 2 pπq is an observable. Here L 2 pπq denotes the weighted L 2 space for the scalar product pφ, ψq π " ş R d φpxqψpxqπpxqdx and the corresponding norm is denoted by }φ} π . A standard methodology for calculating, or, rather, estimating the integral in ( 1) is to construct a stochastic process tXptqu tą0 in R d , e.g. an Itô diffusion process dXptq " f pXptqq dt `σpXptqq dW t

that is ergodic with respect to the measure π. Here W t is a standard m-dimensional Brownian motion and f : R d Ñ R d and σ : R d Ñ R dˆm are assumed smooth and Lipschitz continuous. In particular, π is the unique normalized solution of the stationary Fokker-Plank equation L ˚π " 0, where L ˚is the L 2 pdxq adjoint of the generator Lφ :" f ¨∇φ `1 2 σσ T : ∇ 2 φ of the SDE (2). 5 In what follows we denote by H the L 2 pdxq adjoint of an operator H and by H 7 its L 2 pπq adjoint. Under appropriate assumptions on the drift and diffusion coefficients, we can prove a strong law of large numbers and a central limit theorem as T Ñ 8,

π T pφq :" 1 T ż T 0 φpXptqq dt Ñ πpφq a.e., X 0 " x, (3)
and we have the following convergence in law

? T pπ T pφq ´πpφqq Ñ N p0, σ 2 φ q, (4)
where σ 2 φ denotes the asymptotic variance of the observable φ, given by the Kipnis-Varadhan formula σ 2 φ " xpφ ´πpφq, p´Lq ´1pφ ´πpφqqy π .

(

) 5
Under the assumption that the generator has a spectral gap in L 2 pπq (see for instance [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF]) we have the following exponential convergence

ˇˇEpφpXptqqq ´πpφq ˇˇď Ce ´λt , (6)
where λ ą 0 is the spectral gap of the generator L.

In this paper, we focus on the overdamped Langevin dynamics for sampling [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF],

dXptq " f pXptqq dt `?2 dW t , (7)
where f pxq :" ´∇V pxq and W t is a standard d-dimensional Brownian motion. The invariant measure of (a) is given by πpdxq " Z ´1e ´V pxq dx, where Z " ş R d e ´V pxq dx is the normalization constant and V : R d Ñ R is a smooth confining potential. A question that has attracted considerable attention in recent years is the construction of modified Langevin dynamics that have better sampling properties in comparison to the standard overdamped Langevin dynamics [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF]. Several modifications of the Langevin dynamics [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF] that can be used in order to sample from π are presented in [START_REF] Duncan | Using Perturbed Underdamped Langevin Dynamics to Efficiently Sample from Probability Distributions[END_REF]Sec 2.2]. A well known technique that was first introduced in [START_REF] Hwang | Accelerating Gaussian diffusions[END_REF][START_REF] Hwang | Accelerating diffusions[END_REF] and analyzed in a series of recent papers, e.g. [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF][START_REF] Rey-Bellet | Improving the convergence of reversible samplers[END_REF][START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Duncan | Variance Reduction Using Nonreversible Langevin Samplers[END_REF] for improving the performance of the Langevin sampler [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF], is to introduce in (2) a divergence-free (with respect to the target distribution) drift perturbation g : R d Ñ R d , dXptq " pf pXptqq `gpXptqqqdt `?2 dW t ,

such that divpgπq " 0. (9) We will refer to [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF] as the divergence-free condition. This condition ensures that the SDE (8) has the same invariant measure π as [START_REF] Duncan | Variance Reduction Using Nonreversible Langevin Samplers[END_REF]. We remark that there are infinitely many vector fields g that satisfy [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF]. A complete characterization of all vector fields that satisfy this condition can be found in [START_REF] Hwang | Accelerating diffusions[END_REF]Prop. 2.2].

It is by now a standard, and not difficult to prove, result that nonreversible dynamics exhibits better properties as a sampling scheme, in the sense that the nonreversible perturbation accelerates convergence to equilibrium and reduces the asymptotic variance. The generator of the nonreversible dynamics ( 8) is given by L D " L `A, (10) where L is the generator of (2) A is defined by Aφ " ∇φ ¨g (in the calculations below we will use the notation Aφ " φ 1 g). The drawback of the nonreversible Langevin sampler ( 8) is that, since the generator of the dynamics is a nonselfadjoint operator, a transient, oscillatory phase is introduced. This transient behaviour can be addressed, in principle, by the use of an appropriate splitting numerical scheme [START_REF] Duncan | Nonreversible Langevin samplers: Splitting schemes, analysis and implementation[END_REF].

In this paper, we introduce and analyze an alternative way for perturbing the overdamped reversible Langevin dynamics that is reversible and enjoys all the advantages of the nonreversible sampler (8), while not suffering from the drawback of its oscillatory transient dynamics. The new dynamics is given by the Stratonovitch perturbation

dXptq " f pXptqq dt `gpXptqq ˝?2 dβ t `?2 dW t , (11)
where g satisfies the divergence-free condition and we assume that β t is a one-dimensional standard Wiener process that is independent of W t . 6 For the Stratonovich-perturbed Langevin dynamics [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF] we have the following result. Theorem 1.1 (Reversibility of the perturbed dynamics) Considered the perturbed dynamics (11), were g satisfies the divergence-free condition [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF]. Then the generator of (11) can be written in the form

L S " L `A2 , (12)
and L S is symmetric in L 2 pπq, i.e. L S " L 7 S . As a consequence of Theorem 1.1, the eigenvalues of L S are real, hence there is no transient behaviour of the dynamics. Remark 1 The proposed modified Langevin sampler (11) can be written in the form of general reversible diffusion process 7 (see [START_REF] Pavliotis | Stochastic processes and applications[END_REF]Ch. 4] for the characterization of diffusion processes that are reversible with respect to a given measure):

dXptq " ´pM ∇V qpXptqq dt `pdivM qpXptqq dt `?2DpXptqq d x W t ,
where M " I `gg T " DD T P R dˆd , D " pI, gq P R dˆpd`1q and x W t " pW t , β t q is a standard d `1 dimensional Brownian motion. Theorem 1.2 (Invariant measure preservation) Under the assumptions of Theorem 1.1, the perturbed dynamics [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF] is ergodic with respect to the measure πpdxq " Z ´1e ´V dx.

Remark 2 We note that Theorems 1.1 and 1.2 remain true for general ergodic SDEs (2) with a Stratonovich perturbation, dXptq " f pXptqq dt `gpXptqq ˝?2 dβ t `σpXptqqdW t , (13) where g is a divergence-free vector field with respect to π and f : R d Ñ R d does not have a gradient structure. This includes, in particular, degenerate diffusions (e.g. when the diffusion matrix σσ T is only positive semidefinite), for example the underdamped Langevin dynamics. Indeed the gradient structure is not used in the proofs. Note however that in the case where the functional form of π is not explicitly known, it can be difficult to compute such a vector field.

6. One can also consider Stratonovich perturbations driven by multidimensional Brownian motions with diffusion functions g j , j " 1, 2, . . . satisfying divpg j πq " 0. A detailed analysis of such perturbed Langevin dynamics will be presented elsewhere.

7. Our results can be extended to cover the case of the preconditioned/Riemannian manifold Markov Chain Monte Carlo Langevin dynamics. The details will be presented elsewhere.

The next theorem shows that, in comparison to the original overdamped Langevin dynamics [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF], the Stratonovich perturbation yields a larger spectral gap and a reduced asymptotic variance. Similarly to the nonreversible deterministic perturbation [START_REF] Lelièvre | Free energy computations[END_REF], this hence leads to an improved reversible sampler for the invariant measure (1), both in terms of speeding up the convergence to equilibrium (6) as well as in terms of reducing the asymptotic variance [START_REF] Hwang | Accelerating Gaussian diffusions[END_REF]. When combined, these results provide us with improved performance when measured in the mean-squared error; see [START_REF] Duncan | Using Perturbed Underdamped Langevin Dynamics to Efficiently Sample from Probability Distributions[END_REF]Sec.2.3].

We recall that, under the assumption that the potential V grows sufficiently fast at infinity, the generator of both the standard Langevin and of the Stratonovich-perturbed dynamics have purely discrete spectrum. Theorem 1.3 (Accelerated convergence and reduced asymptotic variance) Let the assumption of Theorem 1.1 hold and let λ L and λ S denote the gaps of the overdamped Langevin (7) and of the Stratonovich-perturbed dynamics [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF], respectively. Then

λ ď λ S . (14)
Let, furthermore φ P L 2 pπq and denote the corresponding asymptotic variances by σ 2 L pφq and σ 2 S pφq. Then

σ 2 L pφq ě σ 2 S pφq. ( 15
)
Remark 3 When the target distribution is Gaussian, in particular for the two dimensional quadratic potential V pxq " 1 2 px 2 1 `λx 2 2 q with λ ! 1, the standard Langevin dynamics [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF] converges to equilibrium at the very slow rate λ L " λ, and it was shown in [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF] that a perturbation of the form

gpxq " δ θ J∇V pxq, J " `0 1 ´1 0 ˘, (16)
with size δ " λ ´1{2 and θ " 1 yields in the case of a nonreversible perturbation (8) an optimally improved convergence rate λ D " Op1q. For isotropic Gaussians, the optimally reduced asymptotic variance using a nonreversible perturbation can also be calculated [START_REF] Duncan | Variance Reduction Using Nonreversible Langevin Samplers[END_REF]Sec. 4]. Similarly, an improved convergence rate of λ S " Op1q can also be obtained for the reversible perturbation (11) for the same scaling δ " λ ´1{2 and θ " 1{2. Observe that the factor δ θ in (16) yields a perturbation of size Opδq of the Langevin generator L in both perturbed generators L D in (10) and L S in [START_REF] Rey-Bellet | Improving the convergence of reversible samplers[END_REF]. It is important to note that the optimal nonreversible perturbation depends on the optimality criterion used, i.e. on whether our aim is to maximize the rate of convergence to equilibrium or to minimize the asymptotic variance (uniformly over the space of square integrable observables). Contrary to this, the optimal reversible perturbation is the same with respect to these two optimality criteria. This observation will be explored further in a future work together with a complete analysis of optimal Stratonovitch perturbations for Gaussian target distributions.

Proof of the main results

We start by recalling from [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF] that the differential operator A is skew-symmetric in L 2 pπq, i.e. A 7 " ´A. This result follows from an integration by parts and (9). To prove our main results we also use that the original SDE (2) has the the generator Lφ " φ 1 f `∆φ.

(17) Proof of Theorem 1.1 We convert the Stratonovitch SDE ( 11) into an Itô one:

dX " f pXqdt `g1 pXqgpXqdt `gpXq ? 2 dβ t `?2 dW t . ( 18
)
Using the calculation A 2 φ " pφ 1 gq 1 g " φ 1 g 1 g `φ2 pg, gq. we deduce the result ( 12) by applying formula (17) to the SDE ( 18). An immediate consequence of A 7 " ´A is then that pA 2 q 7 " A 2 , i.e. A 2 is L 2 pπq symmetric. As L itself is L 2 pπq symmetric, we have that L S is also L 2 pπq symmetric. l Proof of Theorem 1.2 The L 2 -adjoint satisfies

L S π " L ˚π `A˚p A ˚πq " 0, (19)
where we have used the fact that L ˚π " A ˚π " 0. Hence π is the unique invariant measure of the perturbed dynamics [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF]. l

Proof of Theorem 1.3 We write the generator of the Stratonovich-perturbed dynamics as L S " ´B7 B Á7 A with B " ∇, A " g ¨∇, A 7 " ´A. The quadratic form associated to L S is x´L S φ, φy π " }Bφ} 2 π `}Aφ} 2 π for all φ P H 1 pπq the weighted Sobolev space that is defined in the standard manner. The quadratic form associated to the generator of the reversible Langevin dynamics L " ´B7 B is x´Lφ, φy π " }Bφ} 2 π . Since both L S and L are symmetric operators in L 2 pπq with compact resolvents, the spectral gap of the reversible Langevin dynamics is given by the Rayleigh quotient formula, λ S " min

φPH 1 pπq, ş φπ"0 x´L S φ, φy π }φ} 2 π " min φPH 1 pπq, ş φπ"0 }Bφ} 2 π `}Aφ} 2 π }φ} 2 π ě min φPH 1 pπq, ş φπ"0 }Bφ} 2 π }φ} 2 π " λ L .
To prove the bound on the asymptotic variance, we first write the formula for σ 2 S pφq in the form σ 2

S pφq " xψ S , φy π where ψ S is the solution of the Poisson equation ´LS ψ S " φ, and where without loss of generality we have assumed that ş R d φ π " 0. We also consider ψ L , the solution of the Poisson equation ´Lψ L " φ and using L " L S `A7 A, we obtain π `}Bψ} 2 π ě 0, with ψ :" p´Lq ´1p´A 2 qψ S . l Remark 4 Notice also that the perturbation A 2 is only negative semidefinite. In particular, the null space of the perturbation is (much) larger than that of the generator of the overdamped Langevin dynamics which consists of constants. The amount of improvement in the calculation of the integral in (1) using the long time average depends on the magnitude of the projection of the observable φ on the null space of A 2 . Clearly, if this projection is zero, then the inequality in (15) is strict. The details of these arguments will be presented elsewhere.

σ 2 L pφq

Numerical experiments

In this section, we present some numerical experiments to corroborate our theoretical findings and illustrate the features of the Stratonovitch-perturbed Langevin dynamics [START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF]. Although we are primarily interested in large dimensional problems, we consider for simplicity the following warped Gaussian distribution, as considered in [2, Sec. 5.2], with density πpxq " Z ´1e ´V pxq where V pxq is the two-dimensional potential potential V pxq "

x 2 1 100 `px 2 `bx 2 1 ´100bq 2 , where the parameter b " 0.05 is related to how warped the distribution is. For the purposes of this paper, it is sufficient to consider the family of vector fields gpxq " J∇V pxq, J " ´JT , for all constant skew-symmetric matrices J. In particular, we consider the vector field gpxq defined by (16) and we compare the effect of the nonreversible perturbation with θ " 1 in (8) (Figure 1a) and the new reversible Stratonovitch perturbation with θ " 1{2 in (11) (Figure 1b) for several sizes δ " 1, 64, 256 of the perturbation. We also include for reference the results for the standard overdamped Langevin equation [START_REF] Lelievre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF]. We consider the observable φpxq " x 2 1 `x2

2 and consider the estimator 1 M ř M i"1 φpX piq ptqq » EpφpXptqq. We take the initial condition X 0 " p0, 0q and we plot for M " 10 3 independent realisations X piq ptq, i " 1, . . . , M the error | 1 M ř M i"1 φpX piq ptqq ´πpφq| as a function of time t P r0, 4s. The solution is approximated using the simplest Euler-Maruyama method with very small stepsize ∆t " 10 ´5 (considering the Itô formulation (18)). We observe that although the speed of the convergence EpφpXptqq Ñ πpφq as t Ñ 8 is very slow for the standard overdamped Langevin dynamics (see the nearly horizontal black curve for δ " 0), both perturbations lead to an increase in the speed of the convergence to equilibrium (see the transient phase for small time t) while reducing the asymptotic variance (see the equilibrium phase for large time t ě 2 where the oscillations are only due to Monte-Carlo errors), which corroborates Theorem 1.2 and Theorem 1.3. In addition, the Stratonovitch perturbation yields no oscillatory behavior in contrast to the nonreversible one (see Theorem 1.1). This feature renders the new sampling scheme more amenable to efficient numerical methods. This will be explored further in a future study.

  Stratonovitch-perturbed Langevin dynamics[START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF].

Figure 1 .

 1 Figure 1. Error evolution along time of the average over M " 10 3 trajectories of the nonreversible and the Stratonovitch-perturbed Langevin dynamics for different sizes δ " 0, 16, 128, 256 of the perturbation.

  " xφ, ψ L y π " xp´L S qψ S , ψ L y π " xψ S , p´Lqψ L y π ´xA 2 ψ S , ψ L y π " xψ S , φy π `xA 7 Aψ S , ψ L y π " σ 2 S pφq `xAψ S , Aψ L y π . To prove our claim, it is sufficient to show that xAψ S , Aψ L y π ě 0. We calculate, xAψ S , Aψ L y π " xAψ S , Ap´Lq ´1φy π " xAψ S , Ap´Lq ´1pp´Lq `p´A 2 qqψ S y π " xA 7 Aψ S , pI `p´Lq ´1p´A 2 qqψ S y π " }Aψ S } 2 π `xp´A 2 qψ S , p´Lq ´1p´A 2 qqψ S y π " }Aψ S } 2

We assume that the target probability measure has a density with respect to Lebesgue measure. To simplify the notation, we will denote both the measure and the density by π.

For two matrices A and B we use the notation A : B " tracepA T Bq.

(Gilles Vilmart). 1 . Partially supported by the Swiss National Science Foundation, grant No. 200020 172710. 2 . Partially supported by the EPSRC, grants No. EP/P031587/1 EP/L024926/1 EP/L020564/1. 3 . Partially supported by the Swiss National Science Foundation, grants No. 200021 162404, No. 200020 178752 and No. 200020 144313/1. The authors thank Andrew Duncan and Tony Lelièvre for useful discussions. They acknowledge the kind hospitality of the Erwin Schrödinger International Institute for Mathematics and Physics (ESI), where part of this research was developed under the frame of the Thematic Programme Numerical Analysis of Complex PDE Models in the Sciences. Part of the work was done while GP was visiting the Chair of Computational


Abderrahim Mesbah

Hicham Hammouchi

Aissam Berrahou

Hassan Berbia

Hassan Qjidaa

Mohamed Daoudi

Lip Reading with Hahn Convolutional Neural Networks moments

Keywords: Visual speech recognition, Lipreading, Laryngectomy, Hahn

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

moments, Convolutional Neural Networks

Introduction

The visual speech recognition also known by lip reading is a vital task in communication for people with difficulties to interact with society. It has gained a lot of attention lately due to the need of the application of this process in many domains, especially in the medical field which is our motivation for this work. In case of laryngeal cancer, the persons affected by this disease face the loss of their natural voice after the laryngectomy surgery.

Similarly, vocal cord paralysis is another disease that causes the same problems. By assuming that voice therapy can take long to recover the voice, the lipreading is an inevitable task to reconstruct the communication process and save the interaction with society. To concretize that, the advances in computer vision and image processing made this process possible to decode speech from lips movements and offer to these patients a hope to recuperate their communication functions.

Lipreading is a challenging process for humans especially when the context is absent. It requires special qualities for experts to follow lips movements, tongue articulations and teeth. Another confusing issue is the similarity between phonemes explained by Fisher in 1968 [START_REF] Fisher | Confusions Among Visually Perceived Consonants[END_REF]. Additionally, the differences between each speakers mouth shape, mustache, or the effect of makeup can make the task of lipreading more complicated. To face these issues a robust lipreading system is needed to differentiate all these variations.

In the comparison between human and machine lipreading performance conducted by Hilder et al. in 2009 [2], the experiments showed that machine lipreading has outperformed the human lipreading and therefore an automated lipreading system is indispensable to solve the issue.

Toward building an automated lipreading system, several approaches were proposed and tested on several datasets, especially on AVLetters dataset [START_REF] Matthews | Extraction of Visual Features for Lipreading[END_REF], for example an approach using Active Shape Models (ASM) [START_REF] Cootes | Active Shape Models-Their Training and Application[END_REF] and Active Appearance Models (AAM) [START_REF] Cootes | Active Appearance Models[END_REF] was conducted by Matthews et al. in 2002 [5] to extract features from lips images, and train a model using Hidden Markov model (HMM) classifier, this method obtained 44.6% accuracy. Zhao et al. in 2009 [START_REF] Zhao | Lipreading With Local Spatiotemporal Descriptors[END_REF] proposed a lip-reading method using local spatiotemporal descriptors, in which they represented the isolated phrase sequences by extracting the spatiotemporal local binary patterns (LBP) from mouth region.

The best performance attained was 58.85% using Support Vector Machine (SVM) classifier. A method based on Deep bottleneck features extraction directly from pixels was introduced by Petridis and Pantic, 2016 [START_REF] Petridis | Deep Complementary Bottleneck Features for Visual Speech[END_REF], where the authors trained a model using Long-Short Term Memory (LSTM), this method achieved 58.1% accuracy. Bakry and Elgammal in 2016 [START_REF] Bakry | Manifold-Kernels Comparison in MK-PLS for Visual Speech Recognition[END_REF], conducted a comparison between manifold kernels in Manifold Kernel Partial Least Squares(MKPLS). Their approach consists of using distances such Euclid distance between images and LBP to extract features. Another method proposed by Tian and Weijun [START_REF] Tian | Auxiliary Multimodal LSTM for Audio-visual Speech Recognition and Lipreading[END_REF], in which they introduced an auxiliary multimodal LSTM (am-LSTM) that aims to combine audio-visual data at the same time. It learns from both audio and video modalities and uses a pre-trained VGG-16 model to extract features and PCA to reduce dimensionality. On cross modality protocol, which means the audio and video are used for training and only the video is used for testing, the performance obtained was 88.83%.

As for the OuluVS2 dataset, it was first proposed by Annie et al. in 2015 [START_REF] Anina | OuluVS2: A multi-view audiovisual dataset for non-rigid mouth motion analysis[END_REF] to address the problem of non-rigid mouth motion analysis. The provided baseline performance was 41% accuracy on the frontal view. On the same dataset Joon Son Chung and Andrew Zisserman conducted in [START_REF] Chung | Out of Time: Automated Lip Sync in the Wild[END_REF] Further, Saitoh et al. in [START_REF] Saitoh | Concatenated Frame Image Based CNN for Visual Speech Recognition[END_REF] propose a method called CFI-based CNN to represent the spatiotemporal aspect in the video for visual speech recognition. They evaluated the method on OuluVS2 dataset and the performances achieved on the OuluVS2 digits (frontal view) were 61.7% using Netwotk in Network (NiN) model with data augmentation (DA), and an accuracy of 89.4% with DA and using GoogLeNet model. Additionally, regarding the Lip Reading in the Wild (LRW), it was first generated by Chung and Zisserman [START_REF] Chung | Lip Reading in the Wild[END_REF] in 2016, to recognize spoken words. Authors achieved an accuracy of 61.1% using a multiple towers architecture, which uses a convolutional tower for each frame and they concatenate after the first convolution then apply a pretrained VGG-M model. Chung et. al. in [26] present a method called watch, listen, attend and spell, that aims to convert mouth motion videos to characters. They obtained on LRW an accuracy of 76.2% using both CNN and LSTM to recognize the spoken words. Also, Stafylakis and Tzimiropoulos in [START_REF] Stafylakis | Combining Residual Networks with LSTMs for Lipreading[END_REF] and [START_REF] Stafylakis | Deep word embeddings for visual speech recognition[END_REF], obtained 82.97% and 88.08% respectively, This paper is structured as follows. In Section 2 we will introduce Hahn moments concept and mathematical foundation. In Section 3 we will highlight our proposed architecture HCNN. Our Experiments and Results on the reported datasets are found in Section 4. And finally a conclusion to this paper will be given in Section 5.

Two-dimensional Hahn moments

Hahn moments are a set of orthogonal moment based on the discrete Hahn polynomials defined over the image coordinate space. Their implementation does not involve any numerical approximation. In this section, we will give brief formulation of 2D weighted Hahn moments including polynomials and we will describe their capacity to capture significant features from image with significant reduce of dimensionality.

Hahn polynomials

Hahn polynomials of one variable x, with the order n, defined in the interval [0, N -1] as given in [START_REF] Zhou | Image Analysis by Discrete Orthogonal Hahn Moments[END_REF], respect the following equation:

h n (α, β, N |x) = 3 F 2   -n, n + α + β, -x α + 1, -N 1   (1) with n, x = 0, 1, • • • , N -1
where α and β are free parameters, and 3 F 2 is the generalized hyper-geometric function given by :

3 F 2   a 1 , a 2 , a 3 b 1 , b 2 z   = ∞ k=0 (a 1 ) k (a 2 ) k (a 3 ) k (b 1 ) k (b 2 ) k k! z k (2) N -1 x=0 h n α, β, N x h m α, β, N x ω h (x) = ρ h (n)δ mn (3)
where w h (x) is the weighting function given by

ω h (x) = (α + 1) x (β + 1) N -x (N -x)!x! (4)
while ρ h is the squared-norm expressed by

ρ h (n) = (-1) n n!(β + 1) n (α + β + n + 1) N +1 (-N ) n (2n + α + β + 1)N !(α + 1) n (5)
To assure the numerical stability, the set of the weighted Hahn polynomials is defined as

hn α, β, N x = h n α, β, N x w h (x) ρ h (n) (6)
The set of weighted Hahn polynomials obeys the three term recurrence relation defined as follow

hn α, β, N x =A ρ h (n -1) ρ h (n) hn-1 α, β, N x -B ρ h (n -2) ρ h (n) hn-2 α, β, N x (7) n = 2, 3, • • • , N -1
Where

A = 1 + B -x (2n + α + β + 1)(2n + α + β + 2) (n + α + β + 1)(α + n + 1)(N -n) (8) B = n(n + β)(α + β + n + N + 1)(2n + α + β + 2) (2n + α + β)(α + β + n + 1)(α + n + 1)(N -n) (9)
The initial values for the above recursion can be obtained from

h0 α, β, N x = ω h (x) ρ h (0) (10) h1 α, β, N x = 1 - x(α + β + 2) (α + 1)N ω h (x) ρ h (1) (11)

Weighted Hahn moments

The 2D weighted Hahn moments of order m × n of an image intensity function f (x, y) are defined over the domain [0, M -1] × [0, N -1] as:

H mn = M -1 x=o N -1 y=0 hm α, β, M x hn α, β, N x f (x, y) (12)
where M × N is the size of the image. Due to the orthogonality property of separable Hahn polynomials allows us to reconstruct perfectly the image f (x, y), if all moments are used, by utilizing the following inverse transform:

f (x, y) = M -1 x=o N -1 y=0 hm α, β, M x hn α, β, N x H mn (13)
where 0

≤ M ≤ M , 0 ≤ N ≤ N . The image moment set H mn (0 ≤ m ≤ M -1, 0 ≤ n ≤ N -1)
, can be speedily extracted by using the following matrix form [START_REF] Zhou | Image Analysis by Discrete Orthogonal Hahn Moments[END_REF] where H is an m × n matrix of moments H = H ij , by

H = Q T 1 AQ 2 (14)
0 ≤ i ≤ n, 0 ≤ j ≤ m Q 1 = hm (α, β, M x), with 0 ≤ i ≤ m and 0 ≤ x ≤ M -1, Q 2 = hn (α, β, N y), with 0 ≤ j ≤ n and 0 ≤ y ≤ N -1,
and

A = f (x, y), 0 ≤ x ≤ M -1, 0 ≤ y ≤ N -1.
Similarly the image can be reconstructed by

f = Q 1 HQ T 2 (15)
As illustrated in fig. 1, the image can be reconstructed using eq. 15, and obviously despite the small orders Hahn moments preserve the information in image.

Proposed Architecture HCNN

The novel architecture HCNN as shown in fig. 2 on the study conducted in [START_REF] Mesbah | Fast and efficient computation of three-dimensional Hahn moments[END_REF] we have set these parameters to α = β = 5, so the moments retrieved can encompass the whole image with the potential to apprehend its global features. The architecture is divided on two main phases, the Hahn moments layer and the convolutional neural network.

• Hahn Moments layer: the discrete orthogonal Hahn moments as the first layer of HCNN to calculate the moments of the input image and yields a matrix of moments with a size depending on the moments order value.

Therefore, this layer gives an optimized representation of the image and reduce significantly the dimensionality of processing.

• Convolutional Neural Network: comes to further learn pattern in the data and to provide a robust classification. The CNN takes the moments matrix as input instead of the image and apply the various convolutional filters and optimization functions. The first layers of convolution, activation functions, normalization and max pooling, process the input and learn more complex patterns and structures in the data. While the classification is performed in the fully connected layers in which we apply several operations as the normalization, the activation functions and the dropout.

Experiments & Results

In this section, simulation results are carried out through a set of appropriate experiments in order to evaluate the classification performance of the proposed architecture (HCNN). Three datasets have been selected from the literature : AVLetters, OuluVS2 and BBC LRW.

datasets

• AVLetters: a dataset that contains 780 videos for 10 speakers, every speaker utters the 26 alphabet letters three times, which results in 78 videos for each speaker. The speaker start and finish with a closed mouth, and every video has a variable number of frames. The frames as shown in figure4

, are an example of images in AVLetters dataset for the mouth region with a dimension of 80x60.

• OuluVS2: a dataset that contains 52 speakers uttering 10 digits sequences with a repetition of 3 times each. The digits series are as follow : "1 7 3 5 1 6 2 6 6 7", "4 0 2 9 1 8 5 9 0 4", "1 9 0 7 8 8 0 3 2 8", "4 9 1 2 1 1 8 5 5 1", "8 6 3 5 4 0 2 1 1 2", "2 3 9 0 0 1 6 7 6 4", "5 2 7 1 6 1 3 6 7 0", "9 7

Preprocessing

In order to prepare the datasets, we first proceed to retrieve frames from the given videos in both AVLetters and OuluVS2 digits datasets. The extracted frames appear in various numbers depending on video length, which make their manipulation unfitting. To correct this situation we supplement the last frame in each video at the end to unify the numbers of frames in all videos. We choose to add the last frame from each video because we believe that it will not affect the content of the video, since the last frame is just repetition and does not add any information. For the AVLetters we have set the number of frames to 40 frames in each video, so for each video which the number of frames is under 40 we repeat the last frame until we get 40 frames. In the case of OuluVS2 dataset the number of frames is relatively high (the sequence length ranges from 130 to 220), so we proceed differently.

We process each video separately depending on its sequence length and we repeat the last frame until we get a perfect square number of frames, and to determine that number we use the following manner:

let n f be the number of frames in a video, the number c of images to be concatenated as given as follow :

c = ( √ n f ) 2 by x = min{n ∈ Z | n ≥ x}
A further issue to tackle is the modeling of spatiotemporal aspect in the video.

For this matter we use the method of concatenated frame image (CFI) proposed by Saitoh et al. in [START_REF] Saitoh | Concatenated Frame Image Based CNN for Visual Speech Recognition[END_REF] to represent all frames in one concatenated image. The concatenated image is constructed in a way to conserve the chronological order of appearance of frames in the video. For the AVLetters we arrange the 40 frames as one image with 8 frame in each row as illustrated in fig. 7 and fig. 8. In the case of OuluVS2 digits we constructed a square shape image with the size of √ c frames in rows and √ c in columns.

Consequently, we have augmented the OuluVS2 digits dataset to 7650 images.

As for the Lip Reading in the Wild dataset, we start by retrieving the frames from the given videos, then we use the provided meta data to cut the video and keep only the frames corresponding to the spoken word. finally we crop the region of each speaker's mouth (ROI).

Model training parameters

In this section, we highlight the different parameters we used to train our model. After the preprocessing discussed in data augmentation and preprocessing section above, we split our dataset into train and test. For the AVLetters dataset we followed the split given in [START_REF] Matthews | Extraction of Visual Features for Lipreading[END_REF] by taking the two utterances of each speaker for the training and the remaining utterance for the test, which results in 520 images for training and 260 for test, and we train our model for 5000 epochs using a batch size of 520. The OuluVS2 dataset we used the speaker-independent protocol where the twelve speakers (s06, s08, s09, s15, s26, s30, s34, s43, s44, s49, s51, and s52, 10 males and 2 females) were used for testing and the remaining speakers for the training.

After the split we had 1800 images for test and 5850 for train, used to train the model for 1000 iterations with a batch size of 1070. As for the Lip To train our model we first apply the moments filter on each CFI to generate its moments matrix at different orders, where the choice of orders is based on the concatenated image size, as for the images of AVLetters dataset we have calculated moments up to the 16 th , 32 th , 52 th and 60 th , 64 th and 72 th order. In the OuluVS2 as the generated concatenated images are large we have computed the moments up to the 12 th , 16 th , 32 th , 44 th , 56 th and 60 th . The choice of moments order is based on the performance obtained on the test set. We change the order incrementally and we report the results obtained using several orders to show how the performance changes. In addition, We perform an ablation study by removing the Hahn filter and use the same architecture parameters reported in table 1. For AVLetters we resize the CFI to (150 × 320), and for OuluVS2 we set the CFI size to (225 × 225).

In the case of LRW dataset, we use the same model as above (2D for the image dimensions and the third dimension for the sequence length). As we used moments up to the 22 th order, the input has a size of (sequence length× 22 × 22). We feed then the moments matrix to the CNN with parameters illustrated in table 1, with a slight change in the fully connected layers.

We augment the two layers with 300 and 240 neurons to 2048 and 1000 respectively because of the large number of classes (500 words).

Results and Discussion

The recognition rate of our model in comparison with the previous works on AVLetters dataset are shown in table 4. Our method clearly perform better than the methods compared to, which shows the effectiveness of using to represent a whole sequence. Indeed, using Hahn moments we achieve 20% absolute improvement over CNN. As for the OuluVS2 fixed digits dataset we report in table 3 two works for comparison, SyncNet [START_REF] Chung | Out of Time: Automated Lip Sync in the Wild[END_REF] and CFI-based CNN [START_REF] Saitoh | Concatenated Frame Image Based CNN for Visual Speech Recognition[END_REF], where the first employ both CNN and LSTM for recognition Moreover, the results obtained on LRW dataset are listed in table 6. Our method achieves 56.44% top@1% accuracy with data augmentation, with an improvement of about 9% over HCNN without any data augmentation.

Therefore, the augmentation enhances clearly the performance of the model.

In comparison to other works, We perform worse, however, conducting any comparison is not fair, because our model is shallow against the proposed in [START_REF] Chung | Lip Reading in the Wild[END_REF][START_REF] Chung | Lip Reading Sentences in the Wild[END_REF][START_REF] Stafylakis | Combining Residual Networks with LSTMs for Lipreading[END_REF][START_REF] Stafylakis | Deep word embeddings for visual speech recognition[END_REF]. Also, because the data augmentation they applied is extensive against only flipping in our experiments, which allows the model to learn more patterns in the sequences. architectures, despite the shortfalls of some results. These findings would assist researchers to tackle the problem of lipreading and could be a useful aid for the laryngectomized persons to decode their speech from their lips movements. Nevertheless, we believe our work could be a starting point to apply this method on real-time basis, experiment it on real patients, and also to recognize speech independently of the language.

  a method called SyncNet to synchronize mouth motion and speech in a video. The proposed model is a mixture of LSTM and CNN, where the LSTM model ingests the visual features produced by the CNN image by image until the end of the sequence. The model was 92.8% accurate on OuluVS2 fixed digits.

  using a combination of a spatiotemporal 3D CNN, ResNet and Bidirectional The proposed methods above apply spatiotemporal modeling features and visual features extraction before performing classification, which make them computationally expensive and time consuming. In order to cut with these limitations we propose in this paper a novel architecture called HCNN based on Hahn moments and convolutional neural network (CNN). The new architecture lies on the mixture of Hahn moments with its ability to hold and extract the most useful information in images with effectiveness and less redundancy, and the performance of the convolutional neural networks in learning pattern and image classification. The Hahn moments are used as the first layer of our architecture to extract the moments and feed them to the CNN. To the best of our knowledge this is the first time, moments will be used as a filter in a CNN architecture applied to lipreading. The Hahn moments were chosen over other discrete orthogonal moments like chebyshev and Krawtchouk moments, because Hahn moments cover all the properties of both moments, and because of their great ability to represent image with less redundancy in the amount of information. Furthermore, they can be parameterized to retain the global or local characteristics of the image in the lowest orders. In this work we propose a solution that encompasses several issues. The main contributions of this paper are: 1) further improve the performance of the CNN architecture and customize it for a better features extraction and better patterns learning, 2) reduce significantly the dimensionality of images by integrating the Hahn moments as first layer which leads to decrease the computational cost, 3) present a cost-effective solution to the Lip reading problem.

Figure 1 :

 1 Figure 1: image reconstruction with Hahn moments

  aims to solve the problem of Lip reading by processing and recognizing lips images rapidly and efficiently. It is a combination of the method of discrete orthogonal moments and the convolutional neural network. The HCNN architecture comes to surmounts the high computational costs and the sophisticated hardware resources required by the CNN. Furthermore, HCNN enhance the quality of features extraction and the assimilation of patterns incorporated in the image. Indeed, the Hahn moments as a powerful descriptors to retrieve the most useful information in the image, with the property of covering global, local and both features at the same time with efficiency. This advantage can be achieved by tuning the suitable values of α and β parameters as detailed in the Hahn moments section above. Based

Figure 2 :

 2 Figure 2: HCNN: the proposed architecture for lip reading which takes an input image and compute its correspond Hahn moments, then the returned representation is passed through several convolutions, normalization, max pooling, and fully connected layers

Figure 3 :

 3 Figure 3: HCNN model parameters: Hahn moments until the given order, first convolution (kernel 3x3 and 100 filters), second convolution (kernel 3x3 and 60 filters), first max pooling (pool size 3x3). Third convolution (kernel 3x3 and 40 filters). Second max pooling (pool size 3x3). First fully connected layer (300 neurons), and a second layer (240 neurons), and finally an output layer (26 classes)

4 . 2 .

 42 . The dataset as shown in figure 6 is provided with cropped mouth region and with multiple views, with an angle of rotation of the speaker. In our experiments we use the frontal view, with resizing the extracted images to the size of 50×50. • Lip Reading in the Wild: The BBC Lip Reading in the Wild (LRW) dataset contains 500 unique words with up to 1000 utterances per word spoken by different speakers. The dataset as given provide the train, validation and test sets, as well as the metadata indicating the time where the word appear. It is a very challenging dataset with a high number of classes and with words close in spelling which increase the chances of confusion between words. Data Augmentation and Preprocessing 4.2.1. Data Augmentation In order to conduct a fair comparison with other works in the literature, especially those who worked on OuluVS2 dataset, we performed a data augmentation (DA) by applying rotations with angle degrees [-15, -10, 10, 15], on each frame of each video. As for the Lip Reading in the Wild dataset, we augment by flipping horizontally each frame in the dataset. Consequently, the size of the dataset is doubled, and since the dataset is very large, augmenting it by applying other operations makes the processing and the training expensive.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: AVLatters Frames

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Example of CFI for the AVLetters dataset

  Wild dataset we respect the split provided by default, and we train a model for 105 epochs with a batch size of 300.

  in a speaker independent (SI) manner, while the CFI-based CNN lies on a frames concatenation method for modeling the sequences and uses very deep pre-defined CNN architecture such as GoogleNet, AlexNet and Network in Network (NIN) for the recognition task. It can be clearly seen that our shallow HCNN model outperforms the two related works in terms of classification accuracy and reduces enormously the complexity. Similarly to AVLetters dataset, adding Hahn moments filter achieves over than 50% improvement over CNN. Furthermore, in our experiments with only rotation data augmentation we can achieve better than CFI-based CNN and SyncNet, in which, extensive data augmentation like translation, rotation, flipping and color shift were used.

  vides a powerful solution to overcome the highly computation requirements of CNN and deep learning, and to extract the main and useful characteristics of the image to perform the classification with effectiveness. With a small ar-chitecture such the HCNN, we have demonstrated the effectiveness of HCNN on three datasets, AVLetters, OuluVS2 digits and the LRW against deep

Table 1 :

 1 HCNN Model parameters

Table 5 :

 5 Obtained results on OuluVS2 Digits (frontal view) in comparison with other

	methods
	Method	Accuracy
	CFI-based CNN (GoogLeNet) +DA [13]	89.4%
	SyncNet (CNN+LSTM) +DA [22]	92.8%
	HCNN +DA (SI)	93.72%
	CNN Without Hahn +DA (SI)	42.27%

Table 6 :

 6 Obtained results on BBC LRW words dataset whole 500 classes in comparisonIn this paper we introduced HCNN, a novel architecture based on Hahn moments and Convolutional Neural Networks. The proposed method pro-

	with other methods
	Method	T op@1 Acc T op@5 Acc T op@10 Acc
	VGG-M [25]	61.1%	-	90.4%
	Watch-Attend-Spell [26]	76.2%	-	-
	ResNet-LSTM [28]	82.97%	96.3%	98.3%
	Bi-LSTM & ResNet [27]	88.08%	96.28%	-
	HCNN (without DA)	55.86%	82.93%	89.95%
	HCNN (+ flip DA )	58.02%	84.54%	90.86%
	5. Conclusion

Acknowledgment

The authors would like to thank Prof. Ismail Berrada for providing us with the Nvidia GTX 1080Ti GPU used in our experiments.

In the light of the above, the combination of Hahn moments and CNN proved its effectiveness on a complex problem such as the Lip reading problem. HCNN by the component of Hahn moments extract only the significant information by capturing the global features in the image. Further, we significantly reduce the dimensionality and the complexity of the model compared to using the standard CNN architecture. Furthermore, with a small architecture such the ours we can achieve results better than other works that used large architecture such as GoogleNet in the case of OuluVS2, or large pre-trained models like VGG-M.


Antoine Muller

Charles Pontonnier

Georges Dumont

Motion-based prediction of hands and feet contact efforts during asymmetric handling tasks