John M L Ebos
email: john.ebos@roswellpark.org

Michalis Mastri

Christina R Lee

Amanda Tracz

John M Hudson

Kristopher Attwood

William R Cruz-Munoz

Christopher Jedeszko

Peter Burns

Robert S Kerbel
email: robert.kerbel@sunnybrook.ca

Georg A Bjarnason

James G Christensen

Understand (biology) Predict and control (clinic) Differential effects of anti-angiogenic therapies between primary tumor and metastases

Keywords: antibodies, neoadjuvant, surgery, tyrosine kinase inhibitors, VEGF Subject Categories Cancer, Vascular Biology & Angiogenesis

Thousands of cancer patients are currently in clinical trials evaluating antiangiogenic therapy in the neoadjuvant setting, which is the treatment of localized primary tumors prior to surgical intervention. The rationale is that shrinking a tumor will improve surgical outcomes and minimize growth of occult micrometastatic diseasethus delaying post-surgical recurrence and improving survival. But approved VEGF pathway inhibitors have not been tested in clinically relevant neoadjuvant models that compare pre-and post-surgical treatment effects. Using mouse models of breast, kidney, and melanoma metastasis, we demonstrate that primary tumor responses to neoadjuvant VEGFR TKI treatment do not consistently correlate with improved post-surgical survival, with survival worsened in certain settings. Similar negative effects did not extend to protein-based VEGF pathway inhibitors and could be reversed with altered dose, surgical timing, and treatment duration, or when VEGFR TKIs are combined with metronomic 'anti-metastatic' chemotherapy regimens. These studies represent the first attempt to recapitulate the complex clinical parameters of neoadjuvant therapy in mice and identify a novel tool to compare systemic antiangiogenic treatment effects on localized and disseminated disease.

• Theoretical framework for description of the process Can mathematical models be of help in oncology?

Predict and control (clinic)

• Predict tumor growth

y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t, ◊) ' ae M (t, ◊) dN dV dt = aV dV dt = aV " -j = min - X Ä y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t, ◊
y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t, ◊) ' ae M (t, ◊) dN dV dt = aV dV dt = aV " -j = min - X Ä y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t,
y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t, ◊) ' ae M (t, ◊) dN dt = ae ≠bt N dV dt = aV dV dt = aV " -j = min - X Ä y j i ≠ M (t j i , -) ä p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i M : R + ◊ R p ae R m (t,
= aV " -j = min - X Ä y j i ≠ M (t j i , -) ä 2 p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-). 2 - Y j i = M (t j i , -j ) + Á j i , Á j i ≥ N (0, ‡ j i ) M : R + ◊ R p ae R m (t, -) ' ae M (t, -) M : R + ◊ R p ae R m (t, ◊) ' ae M (t, ◊) dN dt = ae ≠bt N dN dt = aN dN dt = aN ≠ eC(t)N N (t) = 2 t T N 0 = e ln(2)t T N 0 … ® dN dt = ln(2) T N (t) = aN (t) N (t = 0) = N 0 dV dt = aV Å 1 ≠ V K ã ® dV dt = aV Ä 1 ≠ V K ä V (t = 0) = 1 mm 3

Competition

Tumor growth

What are minimal biological processes able to recover the kinetics of (experimental) tumor growth? : Gompertz Since the viable tumor mass, mv, is the product of the total number of viable cells and the average mass of a cell, mc, we have mv~Nvmc leading to

dP dt = aP Ä 1 ≠ P K ä dN dt = a K P 2 dV dt = aV dV dt = aV " -j = min - X Ä y j i ≠ M (t j i , -) ä 2 p(y) = 1 Ô 2fi e ≠(y≠-). 2 1 Ô 2fi e ≠(y≠-).
dV dt = aV 1 ≠ V K dN dt = aN Å 1 ≠ N K ã dV dt = aV ln Å K V ã … dV dt = -e ≠ dV dt = aV ln Ä K V ä … dV dt = -e ≠-t V dN dt = aN ln Å K N ã 8 > > < > > : dV dt = aV ln Ä K V ä dK dt = bV 2/3 K(t = 0) = K 0
BT ~( Bc mc )mvz( Ec mc ) dmv dt ð5Þ
This first-order differential equation, representing conservation of energy, explicitly links properties of tumor cells (Bc, Ec, and mc) with properties of the whole tumor (BT and mv). Consequently, it provides a simple, but powerful, way to integrate important features and results from different areas of cancer research. Solving this equation to determine tumor growth requires knowledge of how tumor metabolic rate, BT , depends on its viable mass, mv, to which we now turn.

Model for tumor vascular system and the prediction of metabolic rate. Tumor metabolic rate, BT , is proportional to the sum of the rates of cellular fermentation and aerobic respiration. For avascular tumors, BT depends on the diffusion rate of nutrients and oxygen from the surrounding environment [18]. For vascular tumors, BT is proportional to the total blood volume flow rate to the tumor, _ Q Q T , consistent with observations that glucose and oxygen consumption rates vary linearly with blood flow rate [19]. The dependence of _ Q Q T on mv and host mass, M, is determined by the structure, dynamics and interaction of the tumor and host vasculatures. Here, we develop a complete

® dV dt = aV Ä 1 ≠ V K ä V (t = 0) = 1 mm 3 ® dV dt = aV V (t = 0) = 1 mm 3 ® dV dt = aV log Ä K V ä V (t = 0) = V 0 Y = Y T + Á, Á ≥ N (0, 1) = ® ‡Y - T , Y T Ø V m ‡V - m , Y T < V m -= 0.84, V m = 83 mm 3 8 > > < > > : dV dt = a 0 V, t AE • dV dt = a 1 , t>• V (t = 0) = V 0 ( dV dt = aV " V (t = 0) = 1 mm 3 dV dt = -V " dV dt = aV Å 1 ≠ Å V K ã " ã dN dt = aN Å 1 ≠ N K ã dV dt = aV ln Å K V ã … dV dt = -e ≠-

Metastasis

• Breast

• 94% of cases are local or regional at diagnosis but 30% will relapse

Pollard, N Eng J Med, 2016

• Estimate the amount of residual distant disease at diagnosis in order to personalize the adjuvant (chemo)-therapy • (dormancy Chambers and Groom, Nat Rev Cancer, 2002)

• tumor-tumor interactions • (self-seeding Norton, Nat Med, 2001) 60 Chapter 2. Metastasis: biological dynamics at the organism sca multiplication within organ parenchyma (Fig. 3). These successful metastatic cells ("seed") have been likened to a decathlon champion who must be proficient in 10 events, rather than just a few (64). However, some steps in this pro-cess incorporate stochastic elements. Overall, metastasis favors the survival and growth of a few subpopulations of cells that preexist within the parent neoplasm. The current data, especially studies focused on isolated tumor cells,

Table 1. Steps in the metastatic process

Step Description growth. A detailed review of all the mathematical modeling studies that addressed the question of metastatic dissemination and growth is beyond the scope of the present document (and ca be found, in part, elsewhere [Scott et al., 2013b, Clare et al., 2000]) but a few publications biological or clinical interest that the curious reader can consult are: [Slack et al., 1969, Liot et al., 1976a, Saidel et al., 1976, Guiguet et al., 1982, Koscielny et al., 1985, Klein and BartoszyÒsk 1991, Kimmel and Flehinger, 1991, Yorke et al., 1993, Retsky et al., 1997, Koscielny and Tubian 1999, Iwata et al., 2000, Michor et al., 2006, Hanin et al., 2006, Barbolosi et al., 2009, Bethg et al., 2012, Haeno et al., 2012, Newton et al., 2012, Newton et al., 2013, Scott et al., 2014, Sco et al., 2013a, Mehrara et al., 2013, Coumans et al., 2013, Araujo et al., 2014, Bazhenova et a 2014, Brodbeck et al., 2014, Hanin and Bunimovich-Mendrazitsky, 2014, Hanin et al., 2015, Han and Rose, 2016, Poleszczuk et al., 2016].

The mathematical formalism that we develop below has been employed to address the following b ological questions about metastasis, which have the common characteristic of concerning dynami of the disease at the organism scale.

Biological problem 2. What are the qualitative and quantitative di erences among experiment models of metastasis for di erent cancer types? Is the growth of secondary tumors identical to th Talmadge and Fidler, Cancer Res, 2010 Metastasis: a forgotten major player in modeling • The majority of mathematical modeling efforts in oncology are focused on (primary) tumor growth

• Existing models are based on a statistical, biologically agnostic, prediction of survival

Biologically-based model

Diagnosis

-(V ) = mV - dV p dt = (-p ≠ -p ln(V p )) V p 8 > < > : ˆtfl(t, v) + ˆv ((-≠ -ln(v)) vfl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = µV p (t) " fl(0, v) = 0 8 Molecular data (-omics)
M edical imaging is one of the major factors that have informed medical science and treatment. By assessing the characteristics of human tissue noninvasively, imaging is often used in clinical practice for oncologic diagnosis and treatment guidance 1-3 . A key goal of imaging is 'personalized medicine', where treatment is increasingly tailored on the basis of specific characteristics of the patient and their disease 4 .

Much of the discussion of personalized medicine has focused on molecular characterization using genomic and proteomic technologies. However, as tumours are spatially and temporally heterogeneous, these techniques are limited. They require biopsies or invasive surgeries to extract and analyse what are generally small portions of tumour tissue, which do not allow for a complete characterization of the tumour. Imaging has great potential to guide therapy because it can provide a more comprehensive view of the entire tumour and it can be used on an ongoing basis to monitor the development and progression of the disease or its response to therapy. Further, imaging is noninvasive and is already often repeated during treatment in routine practice, on the contrary of genomics or proteomics, which are still challenging to implement into clinical routine.

The most widely used imaging modality in oncology is X-ray computed tomography (CT), which assesses tissue density. Indeed, CT images of lung cancer tumours exhibit strong contrast reflecting differences in the intensity of a tumour on the image, intratumour texture and tumour shape (Fig. 1a).

However, in clinical practice, tumour response to therapy is only measured using one-or two-dimensional descriptors of tumour size (RECIST and WHO, respectively) 5 . Although a change in tumour size can indicate response to therapy, it often does not predict overall or progression free survival 6,7 . Although some investigations have characterized the appearance of a tumour on CT images, these characteristics are typically described subjectively and qualitatively ('moderate heterogeneity', 'highly spiculated', 'large necrotic core'). However, recent advances in image acquisition, standardization and image analysis allow for objective and precise quantitative imaging descriptors that could potentially be used as noninvasive prognostic or predictive biomarkers.

Radiomics is an emerging field that converts imaging data into a high dimensional mineable feature space using a large number of automatically extracted data-characterization algorithms 8,9 . We hypothesize that these imaging features capture distinct phenotypic differences of tumours and may have prognostic power and thus clinical significance across different diseases. Here we assess the clinical relevance of 440 radiomic features, many of which currently have no known clinical significance, in seven independent cohorts consisting of 1,019 lung cancer and headand-neck cancer patients. Two data sets are used to assess the stability of the features, four data sets to assess the prognostic value of radiomic features on lung cancer patients and head-and-neck cancer patients, and one data set for association Imaging data (Radiomics) Aerts et al., Nat Commun, 2014 letters to nature after their initial diagnosis for an interval of at least 5 years (good prognosis group, mean follow-up of 8.7 years), and 34 patients had developed distant metastases within 5 years (poor prognosis group, mean time to metastases 2.5 years) (Fig. 2a). To identify reliably good and poor prognostic tumours, we used a powerful three-step supervised classi®cation method, similar to those used previously 8,17,18 . In brief, approximately 5,000 genes (signi®cantly regulated in more than 3 tumours out of 78) were selected from the 25,000 genes on the microarray. The correlation coef®cient of the expression for each gene with disease outcome was calculated and 231 genes were found to be signi®cantly associated with disease outcome (correlation coef®cient ,-0.3 or .0.3). In the second step, these 231 genes were rank-ordered on the basis of the magnitude of the correlation coef®cient. Third, the number of genes in the `prognosis classi®er' was optimized by sequentially adding subsets of 5 genes from the top of this rank-ordered list and 2a). To identify reliably good and poor prognostic tumours, we used a powerful three-step supervised classi®cation method, similar to those used previously 8,17,18 . In brief, approximately 5,000 genes (signi®cantly regulated in more than 3 tumours out of 78) were selected from the 25,000 genes on the microarray. The correlation coef®cient of the expression for each gene with disease outcome was calculated and 231 genes were found to be signi®cantly associated with disease outcome (correlation coef®cient ,-0.3 or .0.3). In the second step, these 231 genes were rank-ordered on the basis of the magnitude of the correlation coef®cient. Third, the number of genes in the `prognosis classi®er' was optimized by sequentially adding subsets of 5 genes from the top of this rank-ordered list and PT growth law: g p (V p )=V p (α p -β p ln(V p ))

AL080059 Contig 63649RC LOC51203 Contig 46218RC Contig 38288RC AA555029RC Contig 28552RC FLT1 MMP9 DC13 EXT1 AL137718 PK428 HEC ECT2 GMPS Contig 32185RC UCH37 Contig 35251RC KIAA1067 GNAZ SERF1A OXCT ORC6L L2DTL PRC1 AF052162 COL4A2 KIAA0175 RAB6B Contig 55725RC DCK CENPA SM20 MCM6 AKAP2 Contig 56457RC RFC4 DKFZP564D0462 SLC2A3 MP1 Contig 40831RC Contig 24252RC FLJ11190 Contig 51464RC IGFBP5 IGFBP5 CCNE2 ESM1 Contig 20217RC NMU LOC57110 Contig 63102RC PECI AP2B1 CFFM4 PECI TGFB3 Contig 46223RC Contig 55377RC HSA250839 GSTM3 BBC3 CEGP1 Contig 48328RC WISP1 ALDH4 KIAA1442 Contig 32125RC FGF18 -1 0 1 -1 0
AL080059 Contig 63649RC LOC51203 Contig 46218RC Contig 38288RC AA555029RC Contig 28552RC FLT1 MMP9 DC13 EXT1 AL137718 PK428 HEC ECT2 GMPS Contig 32185RC UCH37 Contig 35251RC KIAA1067 GNAZ SERF1A OXCT ORC6L L2DTL PRC1 AF052162 COL4A2 KIAA0175 RAB6B Contig 55725RC DCK CENPA SM20 MCM6 AKAP2 Contig 56457RC RFC4 DKFZP564D0462 SLC2A3 MP1 Contig 40831RC Contig 24252RC FLJ11190 Contig 51464RC IGFBP5 IGFBP5 CCNE2 ESM1 Contig 20217RC NMU LOC57110 Contig 63102RC PECI AP2B1 CFFM4 PECI TGFB3 Contig 46223RC Contig 55377RC HSA250839 GSTM3 BBC3 CEGP1 Contig 48328RC WISP1 ALDH4 KIAA1442 Contig 32125RC FGF18 -1 0 1 -1 0
Metastases growth law: g(V)=V(α-βln(V))

Pre-surgical Post-surgical fast computation of the metastatic burden

M (t) = Z +OE V 0 vfl(t, v)dv = Z t 0 d (V p (t ≠ s)) V (s)ds
• They are spread by the PT with dissemination rate

d(V p (t)) [day -1 ] dt 8 > < > : ˆtfl(t, v) + ˆv ((-≠ -ln(v)) vfl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = µV p (t) " fl(0, v) = 0 ˆtfl(t, V ) + ˆV (g(t, V )fl(t, V )) = 0 t oe]0, +OE[, V oe]V 0 , +OE[ g(t, V 0 )fl(t, V 0 ) = -(V p (t)) + R +OE V 0 -(V )fl(t, V ) dV t oe]0, +OE[ fl(0, v) = fl 0 V oe]V 0 , +OE[ 8 > < > : ˆtfl(t, V ) + ˆV (g(t, V )fl(t, V )) = 0 t oe]0, +OE[, V oe]V 0 , +OE[ g(t, V 0 )fl(t, V 0 ) = -(V p (t)) t oe]0, +OE[ fl(0, v) = fl 0 V oe]V 0 , +OE[ 8 > < > : ˆtfl(t, v) + ˆv(g(v)fl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = d(V p (t)) fl(0, v) = fl 0 g(V 0 )fl(t, V 0 ) = -(V p (t)) + Z +OE V 0 -(V )fl(t, V ) dV g(V 0 )fl(t, V 0 ) = d(V p (t)) = µV p (t) g(V 0 )fl(t, V 0 ) = d(V p (t)) Ç + Z +OE V 0 d(v)fl(t, v)dv å ˆtfl(t, V ) + ˆV (g(t, V )fl(t, V )) = 0 fl(t, V ) ≥ e ⁄ 0 t (V )
(t, V ) + ˆV (g(t, V )fl(t, V )) = 0 t oe]0, +OE[, V oe]V 0 , +OE[ , V 0 )fl(t, V 0 ) = -(V p (t)) t oe]0, +OE[ , v) = fl 0 V oe]V 0 , +OE[ Y _ ] _ [ ˆtfl(t, v) + ˆv(g(v)fl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = d(V p (t)) fl(0, v) = 0 ˆtfl(t, v) + ˆv(g(v)fl(t, v)) = 0 Y _ ] _ [ ˆtfl(t, v) + ˆv(g(v)fl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = d(V p (t)) fl(0, v) = fl 0 g(V 0 )fl(t, V 0 ) = -(V p (t)) + ⁄ +OE V 0 -(V )fl(t, V ) dV V 0 0 -(V ) = mV - dV p dt = (-p ≠ -p ln(V p )) V p dV p dt = g p (V p ), V p (t = 0) = V i dV p dt = g p (t, V p ), V p (t = 0) = V i dV p dt = g p (t, V p ) Y _ ] _ [ ˆtfl(t, v) + ˆv (g(t, v)fl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = µV p (t) fl(0, v) = 0 g(v) = (-≠ -ln(v)) v Y _ ] _ [ ˆtfl(t, v) + ˆv ((-≠ -ln(v)) vfl(t, v)) = 0 g(V 0 )fl(t, V 0 ) = µV p (t) " fl(0, v) = 0
• Secondary tumors grow in size with rate g(v)

Validation on animal data (Folkman, 2007). The purpose of the present study was to test whether a VEGF RTKI, when administered daily for short periods, can influence the growth of experimental and spontaneous metastasis in mice. The rationale for this experimental design was based on a number of recent clinical and preclinical observations. First, while sorafenib and sunitinib have been approved for the treatment of renal cell carcinoma, as well as hepatocellular carcinoma (sorafenib only) and gastrointestinal stromal tumor (sunitinib only), enduring clinical responses are rare. More-over, single-agent use of such drugs, including bevacizumab, has not led to meaningful beneficial activity in many cases.

Z N vis (t) = Z t≠• vis 0 d(V p (t))1 V p (t)AEV d dt M (t) = µ Z t 0 V p (t)V (t ≠ s)ds i = µ pop + -T x i + ÷ i , ÷ i ≥ N (0, Ê 2 ) ◊ i = ◊ pop + -T x i + ÷ i , ÷ i ≥ N (0, Ê 2 ) ◊ i = ◊ pop + ÷ i , ÷ i ≥ N (0, Ê 2 ) ÷ i ≥ N (0, Ê 2 ) i T i
Second, when such agents are administered on a discontinuous schedule, such as with sunitinib (4 weeks on/2 weeks off), tumor regrowth has sometimes been observed during drug-free break periods (Burstein et al., 2008) or when treatment is discontinued (Johannsen et al., 2008). Third, rapid tumor revascularization has been reported when therapy is stopped in preclinical studies (Mancuso et al., 2006). Finally, we recently reported that a number of changes in proangiogenic plasma proteins observed in patients after sunitinib treatment could be recapitulated in non-tumor-bearing mice, and moreover, in a dosedependent manner (Ebos et al., 2007). Together, these results

Cancer Cell

Brief Antiangiogenic Therapy Increases Metastasis

H Control Sunitinib

N vis (t) = ⁄ +OE V vis fl(t, v)dv = ⁄ t≠• vis 0 d(V p (t))dt M (t) = µ ⁄ t 0 V p (t)V (t ≠ s)ds µ i = µ pop + -T x i + ÷ i , ÷ i ≥ N (0, Ê 2 ) Nicolò et al., in preparation, 2019
Mixed-effects statistical model

•

  ) ' ae M (t, ◊) Theoretical framework for description of the process • Test different hypotheses and reject non-valid onesCan mathematical models be of help in oncology?

  ◊) ' ae M (t, ◊)

Figure 4 |

 4 Figure4| Example of PK/PD simulation to optimize a vinorelbine treatment regimen. The empirical metronomic regimen, incorporating a 50 mg fixed dose of vinorelbine on days 1, 3 and 5 (D1-D3-D5 50 mg) of a 7-day cycle (left panels), can provide substantial clinical benefit to many patients; however, mathematical modelling has helped to identify an alternative dynamic dosing schedule (right panels) of 30 mg, 60 mg and 30 mg on days 1, 2 and 4 (D1-D2-D4 30-60-30),

Figure

  Figure 1. Schematic of t tumor supplied by blood sip represents viable tissue, w Schematic of vascular netw model of tumor and host net

1

  After the initial transforming event, the growth of neoplastic cells is progressive and frequently slow; 2 Vascularization is required for a tumor mass to exceed a 1-to 2-mm diameter (200, 201), and the synthesis and secretion of angiogenesis factors has a critical role in establishing a vascular network within the surrounding host tissue (201); 3 Local invasion of the host stroma by tumor cells can occur by multiple mechanisms, including, but not limited to, thin-walled venules and lymphatic channels, both of which offer little resistance to tumor cell invasion (202); 4 Detachment and embolization of tumor cell aggregates, which may be increased in size via interaction with hematopoietic cells within the circulation; 5 Circulation of these emboli within the vascular; both hematologic and lymphatic; 6 Survival of tumor cells that trafficked through the circulation and arrest in a capillary bed; 7 Extravasation of the tumor embolus, by mechanisms similar to those involved in the initial tissue invasion; 8 Proliferation of the tumor cells within the organ parenchyma resulting in a metastatic focus; 9 Establish vascularization, and defenses against host immune responses; and 10Reinitiate these processes for the development of metastases from metastases.

Figure 3 .

 3 Figure3. The process of cancer metastasis consists of sequential, interlinked, and selective steps with some stochastic elements. The outcome of each step is influenced by the interaction of metastatic cellular subpopulations with homeostatic factors. Each step of the metastatic cascade is potentially rate limiting such that failure of a tumor cell to complete any step effectively impedes that portion of the process. Therefore, the formation of clinically relevant metastases represents the survival and growth of selected subpopulations of cells that preexist in primary tumors.

  lymph node status, histological grade, molecular factors, proliferation indices,… van de Vijver et al., A gene expression signature as a predictor of survival in breast cancer, N Eng J Med, 2002 GENE-EXPRESSION SIGNATURE AS A PREDICTOR OF SURVIVAL IN BREAST CANCER law: d(Vp)= PT growth law: gp Metastases growth law: g Pre-diagnosis Post-diagnosis μ(Vp) γ Metastasis: a forgotten major player in modeling Prediction van de Vijver et al., A gene expression signature as a predictor of survival in breast cancer, N Eng J Med, 2002GENE-EXPRESSION SIGNATURE AS A PREDICTOR OF SURVIVAL IN BREAST CANCER

Figure 1 |

 1 Figure 1 | Extracting radiomics data from images. (a) Tumours are different. Example computed tomography (CT) images of lung cancer patients. CT images with tumour contours left, three-dimensional visualizations right. Please note strong phenotypic differences that can be captured with routine CT imaging, such as intratumour heterogeneity and tumour shape. (b) Strategy for extracting radiomics data from images. (I) Experienced physicians contour the tumour areas on all CT slices. (II) Features are extracted from within the defined tumour contours on the CT images, quantifying tumour intensity, shape, texture and wavelet texture. (III) For the analysis the radiomics features are compared with clinical data and gene-expression data.

Figure 2

 2 Figure 2 Supervised classi®cation on prognosis signatures. a, Use of prognostic reporter prognostic classi®er with optimal accuracy; dashed line, with optimized sensitivity. Above

Figure 2

 2 Figure 2 Supervised classi®cation on prognosis signatures. a, Use of prognostic reporter prognostic classi®er with optimal accuracy; dashed line, with optimized sensitivity. Above

  Median model primary tumor 10th and 90th percentiles model primary tumor Data metastatic burden Median model metastatic burden 10th and 90th percentiles model metastatic burden Nonlinear mixed-effects statistical model for interanimal variability vis -P ,ther ,reb µ, -

2

  Benzekry et al. (Ebos), Cancer Res, ⟹ same growth for PT and mets: α p = α, β p = β Ebos lab Roswell Park Cancer Institute Benzekry et al. (Ebos), Cancer Res, 2016 Introduction Eight inhibitors that block the vascular endothelial growth factor (VEGF) pathway have now been approved as first-or second-line treatment in twelve different late-stage cancer types, thus validating antiangiogenesis as a therapeutic modality in treating established metastatic disease and late-stage glioblastoma (Jayson et al, 2012).Stemming from these approvals, several hundred phase II and III trials were initiated to evaluate VEGF pathway inhibitors in earlier stage disease, that is, neoadjuvant (pre-surgical) and adjuvant (postsurgical) treatment settings(Ebos & Kerbel, 2011). Such 'perioperative' treatments are unique in that they typically have defined treatment durations (unlike in late-stage or advanced disease, where treatments are variable depending on response) and are guided by the hypothesis that drug efficacy in advanced metastatic disease would elicit equal or greater improvements in the earlier stages(Tanvetyanon et al, 2005). These benefits-shown with radiation and chemotherapy (Van Cutsem et al, 2009)-would theoretically include control of localized primary cancers which, in turn, would prevent occult micrometastatic disease and improve progression-free survival (PFS)(Ebos & Kerbel, 2011). However, based on recent clinical and preclinical observations, there is growing concern that VEGF pathway inhibitors may not be effective in this setting(Ebos & Kerbel, 2011). First, there have been five failed phase III adjuvant trials with VEGF pathway inhibitors, including four with the VEGF neutralizing antibody bevacizumab (in combination with chemotherapy or an anti-HER2 antibody) in colorectal carcinoma (CRC) (AVANT and C-08) (de Gramont et al, 2012) and triple-negative and HER2 + breast carcinoma (BEATRICE and BETH, respectively) (Cameron et al, 2013), and one with the VEGF receptor tyrosine kinase inhibitor (RTKI) sorafenib in hepatocellular carcinoma (HCC)(Bruix et al, 2014). Second, growing preclinical evidence suggests that unexpected collateral consequences of angiogenesis inhibition may limit efficacy in preventing growth of micrometastatic lesions(Mountzios et al, 2014). Indeed, we and others have demonstrated that VEGF pathway inhibitors can elicit both tumor-and host-mediated reactions to therapy that can offset (reduce) benefits, or even facilitate, early-stage metastatic disease in certain instances(Ebos et al, 2009; Paez-Ribes et al, 2009). Though these latter results have thus far not been confirmed clinically in patients with advanced metastatic disease when therapy is removed(Miles et al, 2010; Blagoev et al, 2013), they underscore a gap in our current understanding of how antiangiogenic therapy may work in different disease stages. They also raise questions about the translational value of preclinical studies in predicting clinical outcomes. This is of immediate concern as few preclinical studies have tested VEGF pathway inhibitors inINTRODUCTIONMultiple strategies for inhibiting the VEGF pathway have been shown in numerous preclinical studies to hinder tumor growth, and the recent approval of a VEGF-neutralizing antibody (bevacizumab) and VEGF receptor tyrosine kinase inhibitors (RTKIs) (sorafenib and sunitinib) has clinically validated targeting VEGF or its receptors (particularly VEGFR2) as an anticancer treatment

Figure 2 .

 2 Photons/s 10 7

  Fig 2: Simulations of effect or no effect of treatment on metastases

  with size larger than the visible size V vis (= 0.5 cm) τ vis = time to reach V vis • Time to relapse (TTR) defined as the time elapsed from diagnosis to the appearance of a first visible metastasis • Parameter β fixed such that carrying capacity = 10 12 cells

  using the saemix R package (SAEM algorithm)Comets, Lavenu, Lavielle, J Stat Softw, 2017 = P(T i Ø t) = e ≠ = P(T i Ø t) = e ≠ ; -i , µ i 22+ Á i , Á i ≥ N (0, ‡ 2 ) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + ++ + + + + + ++ + + ++ + + + + + + ++ + +++ + + + + Nicolò et al., in preparation, 2019

t V Population fit of tumor growth models

	Exponential	Logistic	Gompertz
	A
	B
	C
	D Vaghi et al. (Benzekry), bioRxiv, 2019

Bayesian estimation for prediction of tumor age

	Vaghi et al. (Benzekry), bioRxiv, 2019

With a priori (Bayesian)

	pertz
	)
	ced	No a priori (MLE)
	pertz
	M)
	ced
	pertz
	)
	pertz
	sian
	nce)
	pertz	A	B
	sian
	nce)
		Gompertz
	ced	(LM)
	pertz
	sian
	nce)
	ced
	pertz	Reduced
	sian	Gompertz (LM)
	nce) e 4. Backward predictions computed with likelihood maximization and with Bayesian infer-
	Three examples of backward predictions of individuals A, B and C computed with likelihood maximization

main cause of death (>90%) from solid cancers Lambert and Weinberg, Cell, 2017 Institut Bergonié, Bordeaux Metastasis (µετά = beyond, στάσιζ = place) Some biological questions of interest to mathematical modeling • Minimal model of metastatic dissemination and colonization able to reproduce the systemic dynamics of a solid cancer disease

• Avoid heavy toxicities for low risk patients •

Lung

• 57% of cases are metastatic • Decide whether to use whole brain radiation therapy or just (stereotactic) surgery • Avoid cognitive impairment of the patient Steeg, Nat. Rev. Cancer, 2016 • Metastases are the • Investigate the relevance of several processes: • (early VS late event Klein, Nat Rev Cancer, 2009)

• (metastases of metastases

Gundem et al., Nature, 2015)

effect of therapy Ebos et

	al. (Kerbel), Cancer Cell, 2009
	• ((pre-)metastatic niche Peinado et al. (Lyden), Nature, 2005)
	• systemic inhibition of angiogenesis O'Reilly et al. (Folkman), Cell, 1990s

• (cancer-immune interactions)

• differential

  Primary tumor V p grows with rate g p [size.day -1 ]

	Mathematical formalism
		10 10
	Volume (cells)	10 4 10 6 10 8
				Diagnosis/Surgery
		10 2
		10 0	0	5	15	20	25	30
					Time (years)

• • Population of metastases represented by a density ρ(t,v) [size -1 ] structured in size v

Acknowledgements

• Preclinical data of ortho-surgical animal models of metastases

• Covariates on !:

• KI67 (p = 3e-04)

• HER2 (p = 0.02)

• CD44 (p = 0.1) Mechanistic model with covariates

r u y r I S 5 t f s T 9 h f 4 V V P H g T x 6 p + w M w 7 4 + a D p V 6 + q q K q X G k 2 O w / A 2 m J i c m v 4 y M / u 1 8 e 3 7 j 5 9 z z f m • Covariates on ":

• EGFR (p = 0.035) c-index = 0.62 estimated with 10-folds cross validation

Predicted TTR Observed TTR

Diagnosis personalization

Volume (cells) Breast cancer patient with primary tumor of 4.32 cm

Nothing visible

Chemotherapy personalization • Number of visible metastases (> 10 8 cel.) 5 years after the end of the treatment Adapt the number of cycles to each patient Barbolosi, Verga, 2011


D Nelson-Gruel

P Lanusse

A Oustaloup

V Pommier

DRAFT: ROBUST CONTROL SYSTEM DESIGN FOR MULTIVARIABLE PLANTSWITH LIGHTLY DAMPED MODES

Keywords: lightly damped modes, MIMO, CRONE, control system, fractional order system, frequency domain 1

A robust controller design is proposed for the active suspension system bench-mark problem. The CRONE control system design used is extended to unstable multivariable plants with lightly damped modes and RHP zeros. Decoupling and stabilizing controller K, is achieved for the open-loop transfer matrix. Fractional order transfer functions are used to define all the components of the diagonal open-loop transfer matrix, β. In defining the fractional open-loop transfer function β 0i some elements of the plants, G0 and its inverse must be considered to achieve the stable controller. Optimisation provides the best fractional open-loop βopt. Finally, frequency domain system identification is used to find controller K=G0 -1 βopt.

INTRODUCTION

The aim of the MIMO (multivariable) CRONE CSD (control system design) is to robustify the closed loop dynamic performance through a robust either damping factor, or a robust resonant peak of control. Fractional differentiation is used to define the optimal open-loop. Three generations of CRONE control have been developed, successively extending the application field. The MIMO approach used here is based on the third generation of CRONE CSD (Oustaloup et Mathieu, 1999). It permits treating uncertain and unstable square nxn MIMO plants with RHP (right half plane) zeros and lightly damped modes (resonant frequencies). Section 2 and 3 briefly outlines the extension of the CRONE CSD methodology. Simulation is given in section 4.

CRONE CONTROL METHODOLOGY AND MIMO SYSTEMS

The aim of the CRONE CSD is to find a diagonal open-loop transfer matrix whose n elements are fractional order transfer functions. It is parametered to satisfy the four following objectives:

-perfect decoupling for the nominal plant, -accuracy specifications at low frequencies, -required nominal stability margins of the closed loops (behaviours around the required cut-off frequencies), -specifications on the n control efforts at high frequencies. (1)

After an optimisation of the diagonal open-loop transfer matrix (1), frequency-domain system identification is carried to obtain the fractional controller. Open-loop transfer functions β0i(s) are used to satisfy the three other objectives.

Definition of the diagonal open-loop transfer function elements

Open-loop transfer function behaviours can be described by using the third generation CRONE CSD methodology presented below. This generation of CRONE CSD uses complex non-integer order integration over a choosen frequency range [ωA,ωB]. The complex fractional order, nf=a+ib allows creating a straight line of any direction in the Nichols chart which is called the generalized template figure (1).

figure 1: Generalized template in the Nichols plane.

The real part of nf determines its phase location at frequency ωcg, that is

2 ) ( π f n e -ℜ
, and the imaginary part determines its direction. The generalized template is described by the limitation in the operational plane Cj of the complex noninteger integrator transfer function:

j f C cg                 = n s s i ω β ) ( 0 , (2) with:
-

ω σ j + = s ∈ Cj, - b a n f i + = ∈ Ci.
This transfer function can be described as based on bandlimited complex non-integer integration:

( ) ( ) a b s s C s i         + + = l h sign ω ω β 1 1 0 ( ) b -q b s s C e sign i l h g /i                             + + ℜ × ω ω 1 1 , (3)
with

                            -         = h cg l cg ω ω ω ω arctan arctan b ch C , (4)
and

2 / 1 2 2 1 1                         +         + = h cg l cg g ω ω ω ω C . ( 5
)
The corner frequencies are placed around the extreme frequencies ωA and ωB such that:

h B cg A l ω ω ω ω ω < < < < . ( 6
)
For stable and minimum phase plant the generalized template is taken account in the open-loop transfer function as follows:

( ) ( ) ( )

s s s s i i i i h 0 l 0 ) ( β β β β = , (7) with ( )
i i i n i s C s l l l l         + = 1 ω β , ( 8
)
the order nli fixes the accuracy of each closed-loop, ( )

i i i i n s C s h h h h         + = 1 ω β , (9)
the order nhi permits the elements of the controller to be proper. This third generation CRONE control open-loop transfer function has been defined using the gain-crossover open-loop frequency; however this definition can be made using the closed loop resonant frequency too. Then ωr replace ωcg.

Decoupling and optimised controller

Let G0 be the nominal plant transfer matrix such that

G0(s)=[ ) (s g ij ]i,j∈N and: [ ] N i i i d n diag diag K G i ∈       = = = 0 0 0 β β , (10)
where: -) (s g ij is a strictly proper transfer function,

-N= { } n ,..., 1 , - i i d n i = 0 β
the element of the ith column and row.

As mentioned above the aim of CRONE control for MIMO plants is to find a decoupling controller for the nominal plant. G0 being not diagonal, the problem is to find a decoupling and stabilizing controller K (Vardulkis , 1987). This controller exists if and only if the following hypotheses are true:

[ ] exist 1 0 1 ) ( : - s G H , (11)
[ ] [ ] 0 ) ( ) ( : 0 0 2 = ∩ + + s G P s G Z H , (12) where [
]

) ( 0 s G Z + and [ ] ) ( 0 s G P +
indicate the positive real part zero and pole sets. The controller K(s) is:

N i i i d n diag G G adj G K ∈ -       = = 0 0 0 1 0 ) ( β , ( 13
) with adj(G0(s))=[G0 ij (s)] T =[G0 ji (s)], G0 ij (s) the cofactor corresponding to element ) (s g ij
, and |G0| corresponding to determinant of G0(s).

Thus: i G G k ji ij 0 0 0 β = N j i ∈ ∀ , . (14)
The nominal sensitivity and the complementary sensitivity transfer function matrices are:

( ) ( ) [ ] ( ) [ ] n i i s S diag s I s ≤ ≤ -= + = 1 0 1 0 β S 0 , ( 15
) | β0i (jω)|dB 0 -π 0 -π/2 ωcg ωA ωB -aπ/2 f(b,a) arg β0i (jω) ( ) ( ) [ ] ( ) ( ) [ ] n i i s T diag s s I s ≤ ≤ - = + = 1 0 0 1 0 β β T 0 , (16) with (
) )) ( 1 ( ) ( 0 0 0 s s s T i i i β β + = , (17) ( ) )) ( 1 ( 1 0 0 s s S i i β + = . ( 18
)
For plants other than the nominal, the closed-loop transfer matrices T(s) and S(s) are no longer diagonal. Each diagonal element ) (

) ( 1 ) ( 1 ) ( ) ( s S s S s T s T s ii ii ii ii ii - = - = β . ( 19
)
For each nominal open-loop β0i(s), many generalized templates can border the same required magnitude-contour of the Nichols chart or the same resonant peak i M 0 p . The optimal one minimizes the robustness cost function:

( )

2 1 min max ∑ = - = n i i i M M J p p (20)
where:

• ( )

) ( sup max max ω ω j p ii T G M i = , (21) • (
)

) ( sup min min ω ω j p ii T G M i = , (22)
while respecting the following set of inequalites for ω∈R and i, j∈N:

( ) ( ) ω ω l j ij ij G T T ≥ inf , ( 23
) ( ) ( ) ω ω u j ij ij G T T ≤ sup , ( 24
) ( ) ( ) ω ω u j ij ij G S S ≤ sup , ( 25
) ( ) ( ) ω ω u j ij ij G KS KS ≤ sup , ( 26
) ( ) ( ) ω ω u j ij ij G SG SG ≤ sup , ( 27
)
where G is the nominal and perturbed plant.

As the uncertainties are taken into account by the least conservative method, a non-linear optimization method must be used to find the optimal values of the independent parameters of the fractional open-loop.

CRONE CONTROL DESIGN FOR UNSTABLE MIMO PLANT WITH LIGHTLY DAMPED MODES

Let the nominal plant transfer matrix G0 be:

                    = ) ( ) ( ) ( ) ( ) ( ) ( 1 1 11 0 s g s g s g s g s g s G nn n ij n       , (28) where:
- ) ( ) ( ) ( 0 s h s g s g ij ij ij ⋅ = , - ) (s h ij
is the transfer including the resonant modes and RHP zeros and poles. The inverse of G0 is:

                    = ) ( ) ( ) ( ) ( ) ( ) ( 1 1 11 s p s p s p s p s p s P nn n ij n       , (29)
where:

- ) ( ) ( ) ( 0 s m s p s p ij ij ij ⋅ = , - ) (s m ij
is the transfer function including the resonant modes and RHP zeros and poles.

Thus the controller is:

) ( ) ( ) ( 0 s s p s k i ij ij β ⋅ = , ( 30
) ) ( ) ( ) ( 0 0 s s m s p i ij ij β ⋅ ⋅ = (31)
Now, (31) implies that some RHP zeros, unstable poles, and lightly damped modes of pij must appear in β0i for the controller to be achievable and stable. The aim of this section is to characterize this open-loop transfer matrix.

Lightly damped modes

Some resonant frequencies must be included in the open-loop transfer function β0i for the controller to be achievable and stable.

The first transfer matrix to consider is the input-disturbance sensitivity, T0K -1 . Using (10):

        + ⋅ ⋅ = = - ) ( ) ( ) ( ) ( ) ( 0 0 0 1 0 s n s d s h s g s d G S K T i i ij i ij (32) this transfer function is resonant-free if N j s h s d s d ij i i ∈ ∀ = ) ( ) ( ) ( . ( 33
)
The denominator of the ith open-loop transfer function must satisfy all the following equations:

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 2 1 s h s d s d s h s d s d s h s d s d in i i i i i i i i = = =  (34)
and therefore:

) ( ) ( ) ( s H s d s d i i i = N i ∈ ∀ , (35)
where, transfer functions Hi(s), have in common some lightly damped modes of the ith row of G0.

The second transfer matrix to consider is input sensitivity KS. Using (10):

        + ⋅ ⋅ = = - ) ( ) ( ) ( ) ( ) ( 0 1 0 0 s n s d s n s m s p T G S K j j j ij ij , (36)
this transfer function is resonant-free if

N i s m s n s n ij j j ∈ ∀ = ) ( ) ( ) ( (37)
The numerator of the ith open-loop transfer function must satisfy all the following equations: ) (

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (
) ( ) ( s M s n s n j j j = N j ∈ ∀ , ( 39
)
where transfer functions Mj(s) have in common some lightly damped modes of the jth column of G0 -1 .

Adding some lightly damped modes on the open-loop transfer functions causes resonant frequencies to appear on sensitivity and complementary sensitivity transfer functions. To attenuate their effect, transfer function Qj(s) is included in β0i(s) around each resonant frequency such that:

1 2 1 2 ) ( ' ' 2 ' 2 + +         + +         = j j j j j s s s s s Q ω ξ ω ω ξ ω , ( 40
)
where:

ωj and ωj ' are frequencies close to the resonant frequency, ξ and ξ ' are damping factors,

Unstable Poles and RHP zeros

0 ) ( = a z η ∀ z ∈ C + .
(41) In our case we need: 0

) ( ≥ ij z k η N j i ∈ ∀ , ∀ z ∈ C + (42)
With ( 18), (43)becomes:

[ ]         ≥ i G G k ji z ij z 0 0 0 ) ( β η η , ( 43
) N j i ∈ ∀ , ∀ z ∈ C + . Consequently: ( ) [ ] ( ) ji z z z G G i 0 0 0 ) ( η η η - ≥ β . ( 44
)
The ith open-loop transfer function must satisfy all the following equations:

( ) [ ] ( ) ( ) [ ] ( ) ( ) [ ] ( ) ni z z z i z z z i z z z G G G G G G i i i 0 0 0 2 0 0 0 1 0 0 0 ) ( ) ( ) ( η η η η η η η η η - ≥ - ≥ - ≥ β β β   (45)
and therefore:

( ) ( ) z G i z z i η η η - ≥ 0 0 ) (β N i ∈ ∀ , ( 46
) with ( ) ( ) ( ) ji z N j i G z 0 min η η ∈ = . (47)
Equation ( 48 Finally, when all the zeros, poles, and lightly damped modes to be include in the ith open-loop transfer function are found:

( ) ( ) ( ) ( ) ( ) s s s s s s s i i i i i i i p h z r l β β β β β β β ) ( ) ( 0 0 = , (48) with:
( ) ( )

i k i z k i i s z C i k z z 0 1 β N i z η β - = Π = , ( 49
) ( ) ( ) i k i p k k i p i i s p s p e C i i k p p 0 1 β N jπ i p η β         - + = Π = Ψ - , ( 50
) ( ) ∏ = = φ β 1 ) ( ) ( ) ( j j i i r s Q s M s H C s i r (51)
where:

-φ is the number of lightly damped modes to integrate.

- k i z is a z such that ( ) 0 0 > i z β η , - k i p is a z such that ( ) 0 0 < i z β η , - ( ) ( ) i i k i z z 0 0 β β η η = if z is a zero of β0i(s), - ( ) ( ) i i k i p z 0 0 β β η η - = if z is a pole of β0i(s), - i z N is the number of RHP zeros of β0i(s), - i p N is the number of RHP poles of β0i(s), - ( ) ∑ = = Ψ i p N β 1 0 k p i k i p i η , - ( ) ∑ = = Ψ i z N β 1 0 k z i k i z i η .
In addition, the controller must be proper and permit the rejection of low amplitude disturbances.

Thus using (30), the controller is proper if:

( ) ( ) ( ) ( ) ) ( ) ( ) ( ) ( s n N D s d i s p s p i ij ij deg deg deg deg + ≥ + , ( 52
)
∀ j ∈ N, with ) ( ) ( s p s p ij ij ij D N s p = ) ( . As a consequence ( ) ( ) ( ) ( ) (
) 53) and at high frequencies

) ( ) ( max ) ( ) ( s p s p N j i i ij ij D N s n s d deg deg deg deg - ≥ - ∈ , (
( ) i h i n s d = ) ( deg (54) and ( ) i z i s n Ψ = ) ( deg , (55)
gives

( ) ( ) (
)

) ( ) ( max s p s p N j z h ij ij i i D N n deg deg - + Ψ ≥ ∈ . ( 56
)
At low frequencies, disturbances are rejected if transfer function SG contains no poles at z=0 or

( ) ( ) 0 ) ( ) ( 0 0 0 ≥ + s G s S η η .
(57) and S(s) is equivalent to β0i(s) or ( ) ( )

i i i z l n s s S Ψ + = ≈ ) ( ) ( 0 0 0 β η η . ( 58
) Also: ( ) ( ) ( ) ) ( ) ( ) ( s g s g ij ij D N s G 0 0 0 0 η η η - = , (59) with:
) (

) ( s g s g ij ij ij D N s g = ) ( . ( 60
)
Therefore:

( ) ( )

) ( ) ( s g s g z l ij ij i i D N n 0 0 η η - ≥ Ψ + , (61)
∀ j ∈ N.

Finally:

( ) ( ) ( ) i ij ij N j i z s g s g l D N n Ψ - - ≥ ∈ ) ( ) ( max 0 0 η η . ( 62
)

APPLICATION

Consider the widely studied 2x2 active-suspension system (Bao-Tung Lin, 2006, Qing-Guo Wang, 2002and Reinelt, 2000), given in figure 2.

      ⋅             Λ + + Λ + + =       ) ( ) ( ) ( ) ( 2 1 1 1 2 1 1 1 2 1 2 2 2 1 s U s U K s B s M s B s B K s B s M s Y s Y , ( 63
) where: The elements to include are:

( ) ( ) 2 1 2 2 1 3 2 1 1 4 2 1 s K M K M s M M B s M M + + + + = Λ ( ) 2 1 1 1 1 2 1 K K s B K K B + + + (
an RHP zero at order ( )

1 0 25 . 1 = = i z β η
for both loops at z=1.25rad/s.This RHP zero is common to each column of G0 -1 ; a lightly damped mode at 0.822rad/s with a damping factor of 0.263 for the first loop, this mode is common to the first column of G0 -1 .

b) Spécifications

For all parametric states, the following must be satisfied:

zero steady-state error for both outputs, short as possible settling time, perfect decoupling for both loops,

first overshoot less than 20%. With this specification some elements of the open-loop transfer matrix can be initialized. So, with the maximum relative degree of G0 -1 (-1) and all elements to integrate in the open-loop transfer matrix, nh must be greater than 3 for the first loop, and greater than 1 for the second loop. A zero steady-state error for both loops causes nl to be equal to one. To attenuate the effect of the damped mode and the steady oscillation (figure 3) transfer ( 40 The final controller is found by frequency domain identification

Results

A simulation of the active suspension system is used to evaluate the controller performance. Inputs of the system are unit step and introduced at t=0s (first input), t=150s (second input); disturbances are 0.1 amplitude step and introduced at t=400s (disturbance on the first loop) and t=600s (disturbance on the second loop). The uncertain spring element K 2 varies between 0.5 and 2. The simulation results are shown in figure 5 for the nominal plant, and figure 6 for nominal and reparametered plants. Overshoot is almost equal to the imposed limit, but the decoupling specification is satisfied for the nominal plant and preserved for reparametered plants. Also the offset and the steady oscillation have been removed.

CONCLUSION

The multivariable CRONE approach is based on third generation CRONE control system design. Each open-loop transfer function β 0i of the matrix β(s) is thus fractional.

Frequency domain system identification and linear algebra are used to determine the final controller. For unstable multivariable plants with lightly damped modes and RHP zeros, some RHP poles, RHP zeros, and lightly damped modes of G0 or G0 -1 can appear on K or on closed loop transfer matrices. The controller would not be achievable and stable. So a method is therefore proposed for treating this kind of system and to determine a final stable and achievable controller. Adding some RHP zeros, RHP poles, and lightly damped modes in open-loop transfer functions before their optimization render the controller achievable and stable. The method used for the active suspension bench-mark problem.

The control obtain is rapid and disturbances rejected.

A crone approach to unstable and non-square multivariable plants with RHP zeros, and lightly damped modes, is being investigated.

Bao-Tung Lin, Wen-Nan Huang, Chin-Cheng

  ) shows that the ith open-loop transfer function β0i(s) must have, for each z, an ( ) zero or a unstable pole of β0i(s).

  64)M1=1 and M2=1.85 are the two masses; K1=0.8 and K2=1.25 are the spring elements of the suspension; and B1=0.8 is the damper. K2 is the uncertain element.

Figure 2 :

 2 Figure 2: System model 4.1. Fractional complex definition of the nominal open-loop transfer matrix

-

  Figure 3: Simulation without feedback

Figure 4 :

 4 Figure 4: Equivalent open-loop transfer function (() nominal, (---) reparametered)

Figure 6 :

 6 Figure 5 : Control performance of the active suspension system

  Teng, «Design of controllers for MIMO active-suspension System Using Combined QFT/ ∞ H », 2006 CACS Automatic Control Conference Oustaloup Alain et Mathieu Benoit, «La commande CRONE», édition Hermes, 1999 Qing-Guo Wang, Yu Zhang, Min-Schen Chiu, «Decoupling internal model control for multivariable systems with multiple time delays», Chemical engineering Sciences, 2002. Reinelt Wolfgang, «Robust control of a Two-Mass-Spring System subject to its Input Constraints», American Control Conference, Chicago, USA 28-30 June 2000 Spencer Jr B.F., S.J. Dyke and H.S. Deoskar, «Benchmark Problems in structural control Part I: Active Mass Driver System», Proceedings of the 1997 ASCE Structures Congress, Portland, Oregon, April 13-16,1997. Tao L., W.D. Zhang, D.Y. Gu, «Decoupling control design processes with Time Delays», Acta Automatica Sinica Vol 31 No.6 2005. Vardulkis A. I. G. -«Internal stabilization and decoupling in linear multivariable systems by unity output feedback compensation», IEEE Trans. on Autom. Control, 32, August, 1987.


Stéphane Loisel
email: stephane.loisel@univ-lyon1.fr.

Pierrick Piette
email: pierrick.piette@gmail.com.

Cheng-Hsien Jason Tsai
email: ctsai@nccu.edu.tw

Applying economic measures to lapse risk management with machine learning approaches

Keywords: JEL Classication: C52, C53, G22 lapse, machine learning, SVM, XGBoost, life insurance

Modeling policyholders lapse behaviors is important to a life insurer since lapses aect pricing, reserving, protability, liquidity, risk management, as well as the solvency of the insurer. Lapse risk is indeed the most signicant life underwriting risk according to European Insurance and Occupational Pensions Authority's Quantitative Impact Study QIS5. In this paper, we introduce two advanced machine learning algorithms for lapse modeling. Then we evaluate the performance of dierent algorithms by means of classical statistical accuracy and protability measure. Moreover, we adopt an innovative point of view on the lapse prediction problem that comes from churn management. We transform the classication problem into a regression question and then perform optimization, which is new for lapse risk management. We apply dierent algorithms to a large real-world insurance dataset. Our results show that XGBoost and SVM outperform CART and logistic regression, especially in terms of the economic validation metric. The optimization after transformation brings out signicant and consistent increases in economic gains.

Introduction

Lapse risk is the most signicant risk associated with life insurance when compared with longevity risk, expenses risk, and catastrophe risk. Policyholders of life insurance may choose to surrender their policies at any time for cash values, or opt to stop paying premiums and leave policies to become invalid eventually. Lapses have signicant impacts on the protability, or even on the solvency, of a life insurer as many studies demonstrate. They may reduce expected prots [START_REF] Hwang | Dierentiating Surrender Propensity from Lapse Propensity across Life Insurance Products[END_REF], cause underwriting expenses unrecovered [START_REF] Tsai | The Distributions of Policy Reserves Considering the Policy-Year Structures of Surrender Rates and Expense Ratios[END_REF][START_REF] Pinquet | Commitment and Lapse Behavior in Long-Term Insurance: A Case Study[END_REF], impair the eectiveness of an insurer's asset-liability management (Kim, 2005c;[START_REF] Eling | Research on lapse in life insurance: what has been done and what needs to be done[END_REF] and bring in liquidity threats as experienced by US life insurers in the late 1980s.

When lapses vary with interest rates as suggested by [START_REF] Dar | Interest Rates, the Emergency Fund Hypothesis and Saving through Endowment Policies: Some Empirical Evidence for the U.K[END_REF], [START_REF] Kuo | An Empirical Study on the Lapse Rate: The Cointegration Approach[END_REF], Kim (2005a), Kim (2005b), and [START_REF] Cox | Annuity Lapse Modeling: Tobit or not Tobit ? Society of Actuaries[END_REF], they become even more detrimental to life insurers [START_REF] Tsai | The Distributions of Policy Reserves Considering the Policy-Year Structures of Surrender Rates and Expense Ratios[END_REF]. Many papers argue that the option to surrender a policy for the cash value might account for a large proportion of the policy value, e.g., [START_REF] Albizzati | Interest Rate Risk Management and Valuation of the Surrender Option in Life Insurance Policies[END_REF], [START_REF] Grosen | Fair valuation of life insurance liabilities: The impact of interest rate guarantees, surrender options, and bonus policies[END_REF], [START_REF] Bacinello | Pricing Guaranteed Life Insurance Participating Policies with Annual Premiums and Surrender Option[END_REF], [START_REF] Bauer | Risk-neutral valuation of participating life insurance contracts[END_REF], [START_REF] Gatzert | Assessing the Risk Potential of Premium Payment Options in Participating Life Insurance Contracts[END_REF], and [START_REF] Consiglio | Pricing the Option to Surrender in Incomplete Markets[END_REF]. The above reasoning and nding may be the reasons why the fth Quantitative Impact Study (QIS5), conducted by the European Insurance and Occupational Pensions Authority (EIOPA) in 2011 regarding the implementation of Solvency II, reports that lapse risk accounts for about 50% of the life underwriting risks.

The signicance of lapse risk draws attentions of scholars to study what causes policyholders to lapse their policies. We may classify the literature into being macro-or micro-oriented. Macrooriented papers (e.g., [START_REF] Dar | Interest Rates, the Emergency Fund Hypothesis and Saving through Endowment Policies: Some Empirical Evidence for the U.K[END_REF][START_REF] Kuo | An Empirical Study on the Lapse Rate: The Cointegration Approach[END_REF]Kim, 2005a;Kim, 2005b;[START_REF] Cox | Annuity Lapse Modeling: Tobit or not Tobit ? Society of Actuaries[END_REF] focus on how lapse rates (the proportion of lapsed policies to the total number of sampled policies within a period of time) are aected by environmental variables such as interest rates, unemployment rates, gross domestic product, and returns in capital markets, as well as by company characteristics like size and organizational form.

Micro-oriented papers secure data from insurers on individual policies to investigate the determinants of the lapse propensities/tendencies. The identied determinants include the characteristics of policyholders and the features of life insurance products/policies (see [START_REF] Renshaw | Statistical analysis of life assurance lapses[END_REF][START_REF] Kagraoka | Modeling Insurance Surrenders by the Negative Binomial Model[END_REF][START_REF] Cerchiara | Generalized Linear Models in Life Insurance: Decrements and Risk Factor under Solvency II[END_REF][START_REF] Milhaud | Surrender triggers in life insurance: what main features aect the surrender behavior in a classical economic context[END_REF][START_REF] Pinquet | Commitment and Lapse Behavior in Long-Term Insurance: A Case Study[END_REF]Eling and Kiesenbauer, 2014, among others) . [START_REF] Eling | Research on lapse in life insurance: what has been done and what needs to be done[END_REF] and [START_REF] Campbell | Modeling of Policyholder Behavior for Life Insurance and Annuity Products: A Survey and Literature Review[END_REF] provide extensive reviews of the literature on lapses 1 .

This paper extends the micro-oriented line of literature in two ways. Firstly, we introduce machine learning algorithms including Extreme Gradient Boosting (XGBoost) and Support Vector Machine (SVM) to lapse behavior modeling. These two advanced algorithms have their merits over other approaches used in the literature such as generalized linear models (i.e., binomial and Poisson models and logistic regression), Classication and Regression Tree (CART) analysis, and the proportional hazards model. Secondly, we adopt economic measures in addition to statistical accuracy in evaluating the performance of dierent algorithms. Such an adoption better demonstrates how dierent algorithms may benet the insurer.

Thirdly, we transform the optimization objective from classication accuracy to economic gains to demonstrate the benet of integrating modeling with prot maximization. Such an integration 1 There are some papers on the subject of modeling early terminations that do not t our macro-micro classication on empirical, explanatory studies. They impose specic assumptions on the transition probabilities to early terminations [START_REF] Buchardt | Cash ows and policyholder behaviour in the semi-Markov life insurance setup[END_REF], the early terminations' intensity [START_REF] Barsotti | Lapse risk in life insurance: Correlation and contagion eects among policyholders' behaviors[END_REF], or the early termination rates [START_REF] Loisel | From deterministic to stochastic surrender risk models: Impact of correlation crises on economic capital[END_REF][START_REF] Milhaud | Exogenous and Endogenous Risk Factors Management to Predict Surrender Behaviours[END_REF]. can increase life insurers' protability, improve insurers' customer management through taking preventive measures to reduce lapses, and retain more of the so-called Contractual Service Margin (CSM) in International Financial Reporting Standard (IFRS) 17. It also links us to the literature on churn management and its impact on the customer lifetime value (e.g., [START_REF] Neslin | Defection Detection: Measuring and Understanding the Predictive Accuracy of Customer Churn Models[END_REF][START_REF] Lemmens | Bagging and Boosting Classication Trees to Predict Churn[END_REF][START_REF] Lemmens | Managing Churn to Maximize Prots[END_REF].

The results from applying dierent algorithms to a large dataset consisting of more than six hundred thousand life insurance policies show that XGBoost and SVM outperform CART and logistic regression with respect to statistic accuracy. The results further show that XGBoost is the most robust across training samples.

The advantages of XGBoost and SVM are more apparent with respect to retention gains. The retention gain takes into account the costs of providing incentives to policyholders to reduce their propensities towards lapses, the benets of retaining policies, and the costs of false alarms. XGBoost and SVM generate much higher retention gains than logistic regression and CART do.

Last but not least, we conrm that economic gains can be further enhanced when the optimization is done on a function linked to the gains rather than on statistic accuracies. The resulted retention gains are 126% of those from applying XGBoost to pursue classication accuracies, and the increase in retention gains remains to be signicant under an alternative policyholder retention scheme. An insurer, therefore, should apply robust machine learning algorithms like XGBoost to its economic objective to achieve optimal lapse management.

The organization of the paper is as follows. Section 2 contains explanations about XGBoost and SVM, followed by brief descriptions on CART and logistic regression. In Section 3 we delineate two performance metrics to be used. One is the commonly seen accuracy, i.e., a statistical validation metric, while the other one is an economic metric considering the expected prots and costs of lapse management. We describe the data obtained from a medium-sized life insurer in Section 4.

Section 5 displays the comparison results across the four algorithms in terms of the statistical and economic metrics. We explain how to integrate algorithms with the prot maximization goal at the beginning of Section 6, and then compare the results from optimizing prot objectives with those from optimization statistic accuracy. Section 7 summarizes and concludes the paper.

Binary classication algorithms

The problem that we want to tackle is detecting whether a policyholder will lapse her/his policy or not, i.e., y i t0, 1u Popular predictive models include logistic regression and CART models. More advanced machine learning models that we introduce in this paper are SVM and XGBoost.

XGBoost

XGBoost is an extension of the gradient boosting introduced by [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF]. The gradient boosting tree is an ensemble method, i.e., multiple weak learners h are combined to become a strong learner F in order to achieve a better predictive performance. The following descriptions are summarized from [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF]).

Given a training sample ty i , x i u N 1 where x i R n and y i t0, 1u, one would like to nd a strong learner F ¦ pxq which minimizes a loss function Ψ py, F pxqq:

F ¦ pxq arg min F pxq E y,x Ψ py, F pxqq $ . (2.1)
The strong learner is an additive expansion of weak learners h ¡ x, tR lm u L 1 , ȳlm © that will be a L-terminal node regression tree in our case:

F M pxq M m0 β m h ¡ x, tR lm u L 1 , ȳlm © M m0 L ļ1 β m ȳlm 1 px R lm q , (2.2)
where tR lm u L 1 and ȳlm are the L-disjoint regions and the corresponding split points determined by the mth regression tree, respectively, and beta m are the expansion coecients. This strong learner is estimated through a stage-wise method that begins with an initial guess F 0 pxq. Then the pseudo-residuals for m 1, 2, . . . , M are computed:

ỹim ¡ δΨ py i , F px i qq δF px i q & F pxqF m¡1 pxq . (2.
3)

The regions tR lm u L 1 are obtained by estimating the mth L-terminal node regression tree on the sample tỹ i , x i u N 1 . The product β m ȳlm γ lm is set to optimize the loss function Ψ:

γ lm arg min γ xi R lm Ψ py i , F m¡1 px i q γq.
(2.4)

At the nal stage, the strong learner is updated,

F m pxq F m¡1 pxq ν.γ lm 1 px R lm q , (2.5)
where ν p0, 1s is a shrinkage parameter that controls how much information is used from the new tree.

The gradient boosting tree method may be summarized as the following algorithm extracted from [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF].

Inspired by previous general works on statistical learning, many extensions to the gradient boosting tree method have been developed. The stochastic gradient boosting technique [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF] is based on the same principle as the bagging technique [START_REF] Breiman | Bagging predictors[END_REF]. It introduces randomness in the observation: given a random permutation π of the integers t1, . . . , N u and Ñ N , the new weak learner tree is estimated on the random subsample 2 ỹπpiqm , x πpiq @ Ñ 1

. Another way to inject randomness that has been popularized by [START_REF] Breiman | Bagging predictors[END_REF] is randomly selecting a subspace of the explanatory variables. More specically, given a random permutation π ¦ of integers t1, . . . , nu and ñ n, the new weak learner tree is estimated on tỹ im , P ¦ pxq i u N 1 in which P ¦ pxq 2

x π ¦ p1q , . . . , x π ¦ pñq @ .

To avoid overtting, some extensions follow the general idea of the ridge regression [START_REF] Hoerl | Ridge Regression: Biased Estimation for Nonorthogonal Problems[END_REF] and lasso regression [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF] and adopt the penalized optimization O py, F pxqq Ψ py, F pxqq ΩpF q.

(2.6) Among all the boosting packages that have been developed, the XGBoost system (Chen and

Guestrin, 2016) has become the most popular due to its exibility and computing performances.

It has also become the most popular machine learning algorithm in data science challenges such as Kaggle for structured data. We list the main parameters that need to be tuned, using the package's terminology and the notation of [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF], as follows.

1. nrounds is the number of trees to grow: M ; 2. eta is the shrinkage parameter: M ν;

3. gamma is the regularization parameter which is used in Ω;

4. max_depth is the number of nodes of a tree: L;

5. min_child_weight is the minimal number of observations in a node and min l,m °N i1 1px i R lm q should be higher than this value; 6. subsample is the relative size of the random subsample used in the case of a stochastic gradient boosting: Ñ {N;

7. colsample_bytree is the relative size of the random subspace of explanatory variables selected at each new tree: ñ{n.

Since we are interested in a binary classication in this paper, we use the logistic loss function:

Ψpy, ŷq N i1 y i ln ¡ 1 e ¡ŷ i © p1 ¡ y i q ln 1 e ¡y i ¨% , (2.7)
and the error function as the metric for cross-validation:

errorpy, ŷq °N i1 1 py i $ round pŷ i qq N , (2.8) where round pŷ i q 5 1 if ŷi ¥ 0.5, 0 if ŷi 0.5.

The tuning method that we adopt consists of two nested cross-validations. We rst perform a grid search on the parameters except nrounds with a 2-folds cross-validation (the grid of values is reported in Appendix 1). Then we determine the best nrounds through a 5-folds cross-validation up to 200 for every possible set of parameters in the grid.

SVM

The theory of SVM was introduced in the 1990's by [START_REF] Boser | A Training Algorithm for Optimal Margin Classiers[END_REF] and [START_REF] Cortes | Support-vector networks[END_REF]. It has become a popular algorithm for classication problems and for churn prediction in particular (e.g., [START_REF] Zhao | Customer Churn Prediction Using Improved One-Class Support Vector Machine. Advanced Data Mining and Applications[END_REF][START_REF] Xia | Model of Customer Churn Prediction on Support Vector Machine[END_REF] Its predictive power is rather good compared to other classication algorithms (e.g., [START_REF] Vafeiadis | A comparison of machine learning techniques for customer churn prediction[END_REF][START_REF] Wainer | Comparison of 14 dierent families of classication algorithms on 115 binary datasets[END_REF].

The SVM algorithm can be described by geometrical terms. The main idea is to nd a hyperplane that separates the observation space into two homogeneous subspaces that is as far apart from each other as possible. This solution is dened as the maximum-margin hyper-plane. To deal with misclassications, a soft margin (i.e., a penalty determined by the user) is imposed upton the SVM.

Another way to deal with classication errors is to project the data to a higher-dimensional space through a kernel function. A more complete geometrical description of SVM can be found in Noble (2006).

In the following, we adopt a formula-based description of the SVM by using the notation of [START_REF] Hsu | A Practical Guide to Support Vector Classication[END_REF]. Given a training sample ty i , x i u N 1 where x i R n and y i t 1, ¡1u, the SVM algorithm is the solution of the following optimization problem:

min ω,b,ξ 1 2 ω t ω C N i1 ξ i , (2.9)
with the constraints

y i ω t φpx i q b ¨¥ 1 ¡ ξ i , ξ i ¥ 0.
(2.10)

The separating hyperplane is determined by the orthogonal vector ω and constant b. The soft margin penalty cost is denoted as C. The data may be projected to a higher dimension space by the function φ, and the underlying kernel function is dened by K px i , x j q φpx i q t φpx j q.

In our case we choose to consider the radial basis function kernel (also called RBF kernel) that is the most commonly used in practice and determined by K px i , x j q exp ¡γ x i ¡ x j 2 ¨, (2.11) with γ ¡ 0 being the kernel parameter.

Then we use the e1071 R package [START_REF] Meyer | Misc Functions of Department of Statistics, Probability, Theory Group (Formely : E1071)[END_REF] to implement the SVM algorithm. To tune the SVM parameters pC, γq, we perform a grid search on a 2-folds cross-validation and adopt the misclassication error function as the validation metric. The grid of values is reported in Appendix 2.

CART

CART was rst introduced by [START_REF] Breiman | Classication and Regression Trees[END_REF]. The underlying idea is straight forward:

dening a class by following a list of decision rules on the explanatory variables. To determine these rules, the data space is iteratively separated by binary split into two disjointed subspaces. At each step or node of this top-down construction, the explanatory variable and the dividing point are chosen to minimize the Gini impurity of the node.

More specically, given a node l of N l observations of response y i t0, 1u with i l, the proportion of observations in the node is dened by p l 1 N l °il y i . Then use an algorithm to partition the parent node into two nodes l L and l R by maximizing

I G plq ¡ rI G pl L q ¡ I G pl R qs , (2.12)
where I G is the Gini impurity of the node and computed by

I G plq N l p l p1 ¡ p l q. (2.13)
This construction is applied up to obtaining a node for every observation point. The tree obtained is thus designated as the saturated model. Although tting the response on the training sample perfectly, it generally leads to low predictive performance when applied to new samples.

Hence the tree needs to be pruned, i.e., the number of nal nodes needs to be reduced to increase its predictive power.

Many criteria can be used to prune the tree, e.g., the minimum number of observations in a nal node. We choose L, the number of terminal nodes, that minimizes the misclassication error: errorpy, ŷq °N i1 1 py i $ ŷi q N .

(2.14)

L is estimated by a 10-folds cross-validation methodology. We use the rpart R package [START_REF] Therneau | Recursive partitioning and regression trees[END_REF] to implement CART.

Logistic Regression

The logistic regression is a special case of the generalized linear models [START_REF] Nelder | Generalized Linear Models[END_REF] obtained with the Bernoulli distribution. The goal is to model the probability of a binary event such as the lapse probability p i of the policyholder i. Given a training sample ty i , x i u N 1 where

x i R n and y i t0, 1u, the regression model is specied as:

ln p i 1 ¡ p i β 0 x t i β.
(2.15)

The parameters pβ 0 , βq R ¢ R n can be estimated by the maximum-likelihood method:

L N ¹ i1 £ e x T i β 1 e x T i β y i ¢ 1 1 e x T i β 1¡y i . (2.16)
When applying the estimated logistic regression model to a classication problem, it doesn't directly lead to labeled responses but to estimated probabilities. To determine the forecasted class, we chose the common threshold of 0.5, i.e., ŷ¦ i 5

1 if ŷi ¥ 0.5, 0 if ŷi 0.5.
(2.17)

Validation metrics

For each policy, the observed lapse y i and the forecasted lapse ŷi are binary variables: py i , ŷi q t0, 1u 2 . The four dierent outputs of a binary classication model are named true positive p1, 1q, true negative p0, 0q, false positive p0, 1q and false negative p1, 0q while the number of each case is usually laid out in the so-called confusion matrix. Denote N pj, kq as the coecients of the confusion matrix in which j t0, 1u stands for the observed lapse indicator and k t0, 1u the predicted lapse indicator. Given a set of response variables ty i , ŷi u N 1 , we estimate N pj, kq as:

N pj, kq N i1
1 py i , ŷi kq .

(3.1)

Statistical metric

Based on the confusion matrix, dierent metrics can be developed. We rst focus on the accuracy metric, the ratio of correctly classied predictions over the total number of predictions:

accuracypy, ŷq N p1, 1q N p0, 0q N 1 ¡ errorpy, ŷq.

(3.2)

Economic metric

Although we adopt mathematic algorithms to predict lapses, the risk is an economic issue after all.

We thus would like to analyze and compare the classication algorithms by an economic metric.

More specically, we will estimate the impacts of dierent classication results on the expected prots from policies, also called customer lifetime values. In order to do so, we plan to adopt an economic model inspired by [START_REF] Neslin | Defection Detection: Measuring and Understanding the Predictive Accuracy of Customer Churn Models[END_REF] and [START_REF] Gupta | Modeling Customer Lifetime Value[END_REF].

Suppose that policy i stays Θ i years in the portfolio pΘ i Nq. The protability ratio at time t can be represented by p i,t and the face amount by F i,t . The lifetime value for policy i is computed as:

CLV i Θ i ţ0 p i,t F i,t p1 d t q t , , (3.3)
where d t is the discount rate.

Assuming a deterministic time horizon T N, we ddene the pT 1q-dimensional real vectors p i , F i , r i and d for protability ratios, face amounts, retention probabilities, and interest rates respectively. Given the four vectors, the customer lifetime value is

CLV i pp i , F i , r i , dq T ţ0 p i,t F i,t r i,t p1 d t q t , , (3.4)
The lapse management strategy is modelled by the oer of an incentive δ i R T 1 to policyholder i who is contacted with a cost c. The incentive is accepted with the probability γ i , and the acceptance will change the vector of the probabilities of staying in the portfolio from r i to r ¦ i R T 1 . We further make the following simplifying assumptions:

1. p i are the same for all policies and denoted as p hereafter;

2. δ i are the same for all contacted policies and denoted as δ hereafter;

3. p i,t , F i,t and d t remain constant across time;

4. r i equals to r lapse or r stay p1, 1, . . . , 1q and r lapse is estimated on the dataset and will be given in Section 5.2;

5. if r i r stay , the incentive is accepted with probability γ i 1 and r ¦ i r stay ; 6. if r i r lapse , the incentive is accepted with probability γ i γ and r ¦ i r stay . 2 Policyholders who reject the oers (probability 1 ¡ γ) will lapse their policies, i.e. r ¦ i r lapse .

The application of a segmentation algorithm to the tested samples produces two confusion matrices: one with respect to number of policies while the other in term of face amount. For the latter matrix, we denote F pj, kq as the coecients of the matrix with regard to face amount, where j stands for the indicator of the policyholder's lapse in real life, k the indicator by the algorithm's prediction, and pj, kq t0, 1u 2 . More specically,

F pj, kq N i1 F i .1 py i , ŷi kq , (3.5)
while N is dened in Equation 3.1.

2 These simplications assume that the protability ratio, the incentive, and the probability to accept the incentive is the same across policyholders, respectively. Upon the availability of data, we may compute an expected protability ratio for each policy. The incentive oered to each policyholder can then be set as a function of the policy's protability. The probability of accepting the oer can also be a function of the incentive, but such a function is dicut to estimate in practice. Face amount may be variable for some products, which increases the diculty in estimating the expected protability ratio. The retention probabilities may change with time, and this calls for a dynamic model of lapse propensities

We dene the reference portfolio value (RPV) as the customer lifetime value of all policies if no customer relationship management about lapses are carried out to be: RP V CLV pp, F p0, 0q F p0, 1q, r stay , dq CLV pp, F p1, 0q F p1, 1q, r lapse , dq .

(3.6) Given a segmentation algorithm, we compute the lapse managed portfolio value (LMPV) by LM P V pδ, γ, cq CLV pp, F p0, 0q, r stay , dq CLV pp, F p1, 0q p1 ¡ γqF p1, 1q, r lapse , dq CLV pp ¡ δ, F p0, 1q γF p1, 1q, r stay , dq ¡ c pNp0, 1q N p1, 1qq .

(3.7)

Then we dene the economic metric of the algorithm as the retention gain:

RGpδ, γ, cq LM P V pδ, γ, cq ¡ RP V,

(3.8)
that can be simplied as RGpδ, γ, cq γ rCLV pp ¡ δ, F p1, 1q, r stay , dq ¡ CLV pp, F p1, 1q, r lapse , dqs CLV pδ, F p0, 1q, r stay , dq ¡ c pNp0, 1q N p1, 1qq . We specify several variables based on the literature and the data provided by the insurer as input to the algorithms of Section 2. Firstly we are able to identify from the data the age, gender, and occupation of an insured at the time when the policy was issued. Female is designated as 1 while male 0 for the dummy variable Gender. Then we designate the dummy variable Occupation as 1 for the occupations that the insurers in Taiwan would undertake extra screening/underwriting. The data also record whether the insured is required to have a physical examination when purchasing life insurance and how many non-life policies (health and long-term care) a person are listed as the insured (since a person may purchase multiple policies).

The data also contain the inception date and face amount of each policy. There are three types of policies. The most popular type is traditional policies like term life, whole life, and endowment.

Investment-linked and interest-adjustable types of products appeared in 2000s. We also able to identify whether a policy is a single-premium one or not. There are three cases with regard to participation. It was not until 2004 that insurers were allowed to sell non-participating policies.

Result with respect to statistical and economic metrics

Our focus is on the predictive performance of dierent algorithms. We thus conduct out-of-sample tests using the following procedure. First, we randomly split the dataset D into 10 subsamples tD 1 , . . . , D 10 u of equal size and then train an algorithm on D k , k t1, . . . , u. The estimated model is subsequently applied to the other subsamples to obtain forecasts ŷ of lapses. In the last step, we compare these predictions with the observed lapses y by the validation metric ρpy, ŷq to measure the predictive performance of the algorithm. This procedure enables us to make sure that every observation is used, at some point of an algorithm, as both training and testing samples. It is similar to the k-fold cross-validation technique in which the training subsample is composed of D ¡D k and the testing subsample is set to D k . We use the k-fold cross-validation to tune parameters in training some of the algorithms.

3 Few policies are also sold by independant agents, brokers that we gather in the same category. 4 Paying premiums by automatic transfers from bank accounts or by recurring payments of credit cards is indierent to policyholders. We thus regard these two automatic/recurring payment methods as one. 5 The exchange rate used in the paper is 30 NTD/1 USD.

5.1

Results with respect to the statistical metric

The mean accuracy computed using the above cross-validation procedure is displayed in the Table 3 and Figure 2 for each binary classication algorithm. As expected, the more sophisticated the model is, the more accurate the predictions will be. XGBoost ranks number one, followed by SVM, CART, and logistic regression (LR). XGBoost surpasses logistic regression by 2.24% on average, which represents a signicant improvement of 12,684 correctly classied policies. Moreover, the smallest standard deviation of accuracy of the XGBoost, 0.03%, indicates that XGBoost is less prone to sample selection. This is visible in the box plot of Figure 2.

Results with respect to the economic metric

To evaluate the algorithms by the economic metric, we rst need to specify the parameters of the cash ows model. Since no data is available for us to estimate these parameters, we have to make assumptions. We had conducted sensitivity analyses and conrmed that the comparison results remain the same in general.

The time horizon T is set to 12 years according to the length of the sampling period. We estimate the retention probability vector r lapse from the dataset and obtain the vector displaid in Table 8. Other parameters are set as follows:

the protability ratio p 0.5%; the discount rate d 2%; the cost to contact a policyholder c 10 USD.

We propose two dierent incentive strategies: an aggressive one and a moderate one. The incentive vectors are described in Table 9 Table 9: Incentive strategies. We further assume that the probabilities of accepting the incentives for a would-lapse policyholder are γ 1 20% and γ 2 10% respectively.

The results from comparing dierent classication algorithms by the economic metric with the aggressive incentive strategy are displayed in Table 10 and Figure 3. The winner looks to be XGBoost: it has the highest retention gain with the smallest standard deviation across subsampling.

Figure 3 further illustrates that XGBoost and SVM lead to similar retention gain compared to logistic regression and CART.

Notice that the dierences across the algorithms are wider in terms of the economic metric than the statistical metric. The accuracies of the models are between 76.64% and 78.88%, which means an improvement ratio of 2.9%. The retention gains, on the other hand, range from 2.7 and 5.2 million USD, indicating an enhancement of 96%. Therefore, choosing a good algorithm is more important in terms of economic reality (dollar amount) than by statistical accuracy. It appears that CART produces the lowest retention gain: $2,680,012. This is mostly because CART has the highest false alarm rate (cf. Table 3c) which means oering the incentive to many policyholders who have no intention to lapse their policies. Furthermore, CART leads to the highest contacting cost since it predicts the highest lapses. The prots are thus reduced. Then we look at algorithms' performances when the incentive strategy is moderate and leads to lower acceptance probabilities. The results are displayed in the Table 11 and the Figure 4. We rst notice XGB and SVM remains to be ranked No. 1 and No. 2, respectively. Next we observe that the improvement ratio of the best algorithm over the worst is smaller but remains to be signicant (56%). Thirdly, retention gains are signicantly lower with the moderate incentive strategy. For instant, XGB achieves a gain of 5.2 million dollars with the aggressive incentive strategy but the gain reduces to 3.3 million dollars when incentives oered to policyholders are moderate. Under our assumptions, the company should rather set the aggressive incentive strategy up to optimize her gains. However, in practice, one would need a more complete sensitivity study on the incentive to be oered and the corresponding acceptance probability to fully optimize the lapse management. 6 Optimization on protability instead of classication It is obvious that insurers would not seek to optimize the classication accuracy but focus on economic gains resulted from the classication algorithms when forming a lapse management strategy.

When our aim is to maximize the protability of the lapse management strategy, binary classications might be unsuitable since they are not designed to meet such a need. [START_REF] Ascarza | In Pursuit of Enhanced Customer Retention Management: Review, Key Issues, and Future Directions[END_REF] emphasize the dierence between the at-risk population (e.g., customers with high churn probabilities) and the targeted population (e.g., customers that the company should focus her retention campaign on in order to optimize her prots) from an economic point of view. Along this line of churn literature, [START_REF] Lemmens | Managing Churn to Maximize Prots[END_REF] modify the usual loss function into a prot-based function to optimize economic gains. They obtain a signicantly increase in the expected prot of a retention campaign. Learning from the churn literature, we transform the above classicaton problem into a regression question in this section.

Methodology

Let the new response variable z R j i represents the retention gain or loss resulting from proposing the incentive j t1, 2u (cf. Section 5.2) to policyholder i. More specically, we dene z R j i as z R j i 5 ¡CLV pδ j , F i , r stay , dq ¡ c if y i 0, γ j . rCLV pp ¡ δ j , F i , r stay , dq ¡ CLV pp, F i , r lapse , dqs ¡ c if y i 1.

Ψ z R j , ẑR j ¨ 1 N N i1 z R j i ¡ ẑR j i % 2 , (6.2)
and as the metric for cross-validation.

In the last step, lapse ŷi is forecasted if the estimated gain is positive:

ŷi 5 1 if ẑR j ¥ 0, 0 if ẑR j 0, (6.3)
By this way we can apply the same metrics described in previous sections. Here ŷi is better to be understood as the estimation of the protability about oering an incentive to the policyholder i rather than the forecast on the policyholder's lapse.

The two new classications are denoted as XGB_R1 and XGB_R2, respectively, for applying XGBoost to z R 1 and z R 2 . The tuning method that we apply to estimating the parameters is described in Appendix 3.

Results

Table 12 and Figure 5 display the prediction accuracies. Table 12 shows that XGB_R1 and XGB_R2 produce relatively low mean accuracy of respectively 76.7% and 75.7% While XGB_R2

is clearly the worst model in term of accuracy, XGB_R1 generates similar results to the logistic regression which is the worst binary classication model regarding the accuracy measure. These seemingly unsatised results are understandable since both XGB_R1 and XGB_R2 are not designed to predict whether a policy would be lapsed or not. What they aim for are economic gains. and 99,432) and produce the most false-sense-of-security (resp. 113,948 and 119,405). However, we will see very soon that XGB_R1 and XGB_R2 stand out when we switch focus to retention gain.

Table 15 and Figure 6 show that XGB_R1 generates a signicantly larger average retention gain with the aggressive incentive strategy ($6,586,357) than other algorithms as well as a signicantly lower standard deviation ($53,460). The increase in retention gain is 26% (1.3 million USD) higher than that generated by XGB (the second-best algorithm) and 146% (3.9 million USD) better than that produced by CART. Looking back to Table 13, we see that XGB_R1 leads to reduce the number of false alarms (18,204) in optimizing the retention gain, even if this also reduces the We expect that the benet of integrating the algorithm with the goal is robust across incentive strategies. This is conrmed by the results in Table 16 and Figure 7. XGB_R2 generates retention gain of 3.9 million dollars that is nearly 600 thousand dollars more than that achieved by the second place XGB. The increase in retention gains is 18%. The increases with respect to the commonly seen LR and CART reach 47% and 85%. The results in this section demonstrate the benet of having a specic objective. If senior managers of an insurer are able to specify an objective to be optimized (e.g., maximizing retention gain), the sta should apply an advanced algorithm like XBG directly to such an objective to achieve the optimum. The enhanced gain relative to the case having no specic objective other than classication accuracy can be substantial.

Conclusions

Lapse risk is the most signicant risk associated with life insurance. Lapses may cause losses, reduce expected prots, lead to stringent liquidity, result in mis-pricing, impair the risk management, or even pose solvency threats. Employing a good algorithm to model policyholder lapse behavior is therefore valuable. In this study, we adopt innovative viewpoints on lapse management in addition to introducing machine learning algorithms to lapse prediction. Applying XGBoost and SVM to predicting whether a policyholder will lapse her/his policy is new to the literature. Secondly, we adopt not only a statistical metric in evaluating algorithms' prediction performance but also an economic metric based on customer lifetime value and retention gains.

The goal of classication accuracy has no direct link to the insurer's costs and prots. It thus might lead to a biased strategy [START_REF] Powers | Evaluation: from Precision, Recall and F-measure to ROC, Informedness, Markedness and Correlation[END_REF]. Following the churn literature, we dene a specic validation metric based on the economic gains. This constitutes our third contribution: we are the rst to set up a prot-based loss function so that we may directly optimize the economic gains.

More specically, we change the usual statistical idea of classication to a gain regression in which prots are to be maximized.

The two machine learning algorithms, XGBoost and SVM, perform a little bit better than classic CART and logistic regression in terms of statistical accuracy on a large dataset consisting of more than six hundred thousand life insurance policies with information on policy terms and policyholders' characteristics. XGBoost has another advantage over other algorithms: it is less dependent upon the choice of training samples.

The advantages of XGBoost and SVM are more apparent with respect to retention gains. The retention gains incorporate the costs of providing incentives to policyholders to reduce lapse propensities and the benets of retaining policies. XGBoost and SVM generate much higher retention gains than logistic regression and CART do. For instance, XGBoost produces 1.2 to 2.6 million dollars more economic gains than CART.

In the last section, we demonstrate that the economic gains can be further enhanced when the optimization is done on a function linked to economic gains rather than on statistic accuracies. The results show that the retention gains with an aggressive incentive strategy resulted from XGB_R1 is 126% of those from applying XGBoost to pursue classication accuracies, in particular by reducing the false alarm rates. An insurer should therefore apply advanced machine learning algorithms like colsample_bytree : 0.8.

Then, we determine the best nrounds through a 5-folds cross-validation with this parameter tested up to 1,000.

Figure 1 :

 1 Figure 1: Pseudocode of the Gradient Tree Boosting algorithm[START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF].

  from a medium-size life insurance company in Taiwan that had total assets over 15 billion US dollars at the end of 2013. The data contain 629,331 life insurance policies sold during the period from 1998 to 2013. The data-providing insurer tracked changes in the statuses of policies including death and lapse. The last tracking date is 31/08/2013. 243,152 policies out of all samples were lapsed, and 5,486 insureds died during the sampling period.

Figure 2 :

 2 Figure 2: Box plot of statistic accuracies.

Figure 3 :

 3 Figure 3: Box plot of retention gains with the aggressive strategy.

Figure 4 :

 4 Figure 4: Box plot of retention gains with the moderate strategy.

Figure 5 :

 5 Figure 5: Box plot of statistic accuracies.

Figure 6 :

 6 Figure 6: Box plot of retention gains with the aggressive strategy.

Figure 7 :

 7 Figure 7: Box plot of retention gains with the moderate strategy.

Table 1 :

 1 Descriptive statistics of categorical explanatory variables.

	Variables	Category	Percentage
	Gender
		Female	48
		Male		52
	Occupation
		Tier one	80.5
		Requiring extra screening	19.5
	Physical Examination
		Exempted	96.4
		Required	3.6
	Distribution Channel
		TA		93.9
		BK		3.4
		DM		2.4
		Others	7	0.3
	Premium Payment
		Single premium	3.1
		Non single premium	96.9
	Premium Paying Method
		Insurer	18.8
		B&C		70.8
		P&C		10.4
	Participation
		Non-participating	37.2
		Participating	16.2
		Mandatory participating	46.6
	Product Type
		Interest-adjustable	1.7
		Investment-linked	1.2
		Traditional	97.1
	Currency Domination
		NTD		88.1
		Others		11.9

Table 2 :

 2 Descriptive statistics of continuous explanatory variables.

		Mean	Medium	St. Dev.	Minimum	Maximum
	Age	28.3	27	16.8	0	80
	# of non-life policies	1.2	0	2	0	33
	Inception date	06/06/2005	21/04/2005	4.8 (years)	01/01/1998	31/07/2013
	Face Amounts (in USD)	17,165	10,000	28,050	333	2,000,000

Table 3 :

 3 Cross-Validated Statistic Accuracies.

		LR	CART	SVM	XGB
	Mean Accuracy	76.64%	77.15%	77.82%	78.8%
	Standard Deviation	0.07%	0.10%	0.08%	0.03%

Table 4 :

 4 Average confustion matrix of XGB.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	309,111 81,177	38,450 137,660

Table 5 :

 5 Average confustion matrix of SVM.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	310,258 88,339	37,303 130,498

Table 6 :

 6 Average confustion matrix of CART.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	296,320 78,209	51,241 140,628

Table 7 :

 7 Average confustion matrix of LR.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	304,025 88,775	43,537 130,062

Table 8 :

 8 Estimated retention probability r lapse .

	Year	0	1	2	3	4	5	6	7	8	9	10	11	12
	Retention probability	0.96	0.87	0.67	0.37	0.27	0.21	0.15	0.12	0.10	0.08	0.06	0.05	0.04

Table 10 :

 10 Cross-validated retention gains with the aggressive strategy.

		LR	CART	SVM	XGB
	Mean Retention Gain	4,046,602	2,680,012	5,028,737	5,243,913
	Standard Deviation	133,993	209,220	139,102	115,415

Table 11 :

 11 Cross-validated retention gains with the moderate strategy. The drawbacks of XGB and SVM relative to CART and LR that we may think of are not

		LR	CART	SVM	XGB
	Mean Retention Gain	2,618,396	2,085,599	3,113,900	3,261,029
	Standard Deviation	63,693	85,184	54,169	45,928
	In summary, XGB and SVM consistently perform better than CART and LR no matter which
	performance index, statistical accuracy or retention gains with alternative incentive strategies, is
	used.

Table 12 :

 12 Cross-Validated Statistic Accuracies.

		LR	CART	SVM	XGB	XGB_R1	XGB_R2
	Mean Accuracy	76.64%	77.15%	77.82%	78.8%	76.67%	75.71%
	Standard Deviation	0.07%	0.10%	0.08%	0.03%	0.07%	0.06%
	The numbers in Table 13 and 14 tell us more about why XGB_R1 and XGB_R2 performs badly
	in statistical accuracy. They result in the smallest correct identications on lapses (resp. 104,889

Table 13 :

 13 Average confustion matrix of XGB_R1.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	329,357 113,948	18,204 104,889

Table 14 :

 14 Average confustion matrix of XGB_R1.

			Predicted
			Stay	Lapse
	Actual	Stay Lapse	329,413 119,405	18,149 99,432
	correct detection (104,889). The good results of XGB_R1 in achieving retention gain demonstrate
	the benet of integrating the algorithm with the goal to be achieved. The objective function for
	XGB_R1 to minimize, Equation 6.2, is about predicting retention gains. XGB_R1 therefore would
	naturally perform the best when compared with other algorithms optimizing other objectives (such
	as classication accuracies).

Table 15 :

 15 Cross-validated retention gains with the aggressive strategy.

		LR	CART	SVM	XGB	XGB_R1
	Mean Retention Gain	4,046,602	2,680,012	5,028,737	5,243,913	6,586,357
	Standard Deviation	133,993	209,220	139,102	115,415	53,460

Table 16 :

 16 Cross-validated retention gains with the moderate strategy.

		LR	CART	SVM	XGB	XGB_R2
	Mean Retention Gain	2,618,396	2,085,599	3,113,900	3,261,029	3,852,782
	Standard Deviation	63,693	85,184	54,169	45,928	39,163

This policy is a whole life insurance with a one-year old insured and the death benet of ten thousand NTD (a little over three hundred USD). There are other small policies with death benets smaller than three thousand USD. These policies constitute less than one percent of our samples.

The author acknowledges support from research chair DAMI (Data Analytics and Models for Insurance) sponsored by BNP Paribas Cardif. The author is grateful to the Ministry of Science and Technology of Taiwan (project number MOST 105-2410-H-004 -070 -MY3) and Risk and Insurance Research Center of National Chengchi University for the nancial supports.

Appendix 1 XGBoost Tuning -Binary Classication

The values of the parameters tested in the grid search for the tuning of XGBoost are as follows: eta : 0.005, 0.1, 0.15; gamma : 0, 5, 10; max_depth : 10, 15, 20, 25, 30; min_child_weight : 15, 20, 25; subsample : 1; colsample_bytree : 0.4, 0.5, 0.6.

The values of the grid search are chosen by a previous sensitivity study in which we apply the same methodology on a subsample of the whole database but with a coarser grid. Then we focus on a ner grid to obtain better results within a reasonable time period. In addition, the fact that we only test subsample with the value of 1 means that we do not adopt the stochastic gradient boosting of [START_REF] Friedman | Greedy Function Approximation: A Gradient Boosting Machine[END_REF].

Appendix 2 SVM Tunning

The values of the parameters tested in the grid search for the tuning of SVM are as follows:

Cost : 0.5, 1, 2, 5, 10; gamma : 0.25, 0.5, 0.75, 1, 1.25.

Similar to the previous section, the values of the grid search are chosen by a previous sensitivity study in which we apply the same methodology on a subsample of the whole database but with a coarser grid. Then we focus on a ner grid to obtain better results. This is necessary so that the computing can be done within a reasonable time period.

Appendix 3 XGBoost Tuning -Protability

We adopt the values of most parameters generated by a previous sensitivity study as: eta : 0.005; gamma : 1; max_depth : 15; min_child_weight : 15; subsample : 0.7;