Jean Clairambault
email: jean.clairambault@inria.fr

Mathematical modelling of cancer growth and drug treatments: taking into account cell population heterogeneity and plasticity

Mathematical models of cancer growth and evolution of cancer cell characteristics, aka traits or phenotypes, together with optimisation and optimal control methods to contain them, in the framework of adaptive dynamics of cell populations, are presented. They take into account the heterogeneity of cancer cell populations, i.e., their biological variability, and their intrinsic plasticity, i.e., their nongenetic instability that allows them to quickly adapt to changing environments. The presented vision of the cancer disease, which is specific to multicellular organisms, relies on a relatively novel vision, consistent with a billion-year evolutionary perspective. Based on recent contributions from philosophy of cancer, these mathematical models aim at designing theoretical therapeutic strategies to simultaneously contain tumour progression and limit adverse events of drugs to healthy cell populations.

I. INTRODUCTION

"What is cancer?" is a question that must somehow be answered with clear ideas before proposing relevant mathematical models of cancer growth and phenotype evolution, together with methods to optimise its treatments. Recent contributions from philosophy of biology, immunology, and archeology of genomes, are reported here, stressing the importance of a default of cohesion within and between tissues of a same organism -as cancer is a disease of multicellular organisms only -, that is manifested as an impairment of organismic control on both proliferation and differentiation.

Different mathematical models of tumour growth involving cancer cell population heterogeneity and plasticity, and drug treatments of cancer are presented in the present paper. The ideas developed here are not a review of the vast domain of studies dedicated to modelling heterogeneous cell populations. Devoted to advocating and illustrating by a few chosen models the interest of continuous phenotypestructured equations in cancer and its treatment, they are intended to be presented as part of a tutorial session in ECC 2023 on "Modelling and Treatment of Cancer". They concern deterministic models, with a focus set on adaptive dynamics of cell populations, healthy and cancer, and the effects on them of anticancer drugs. Constraints linked to avoiding the main pitfalls of therapeutics in oncology, namely therapy resistance in cancer cells and adverse toxic events in healthy cells, together with an objective function that is most often limited to maximising tumour cell kill, are considered, paving the way for the use of optimisation methods. Optimised clinical strategies for the continuous delivery of anticancer drugs may thus be designed by resorting to optimal control methods. Such optimised strategies should thwart some known pathological strategies used by cancer cell populations to dodge deadly insults induced by drugs, such as the development of drug-resistant cell phenotypes and the so-called 'bet hedging' by which tumours may simultaneously face different modes of insults in abruptly changing environments [START_REF] Brutovsky | Structure of Intratumor Heterogeneity: Is Cancer Hedging Its Bets?[END_REF], [START_REF] Gravenmier | Adaptation to Stochastic Temporal Variations in Intratumoral Blood Flow: The Warburg Effect as a Bet Hedging Strategy[END_REF], [START_REF] Philippi | Hedging one's evolutionary bets, revisited[END_REF].

From a modelling point of view, the complete reversibility of the cell phenotypes under study when the environment changes is represented by the reversible behaviour of the solutions of the equations in terms of cell-structuring traits (aka phenotypes). It is consistent with the fact that we study in an isogenic context cell differentiations and dedifferentiations, i.e., reversal of differentiations within a given cell lineage toward a less mature phenotype [START_REF] Shen | Cell plasticity in cancer cell populations (review) [version 1; peer review: 2 approved[END_REF], that rely on epimutations, i.e., epigenetic changes in the expression of genes, and not on gene mutations [START_REF] Huang | Genetic and non-genetic instability in tumor progression: link between the fitness landscape and the epigenetic landscape of cancer cells[END_REF]. With or without interactions between cell populations, the equations are of the nonlocal Lotka-Volterra type, structured by adaptive continuous cell phenotypes, with optional diffusion and advection terms (two examples are given below). Adaptation, i.e., phenotype changes under environmental pressure, e.g., due to competing cell populations, drugs or other therapeutic means, is thus considered at a cell population status at which the situation of a cancer cell population is reversible (i.e., not after mutations), in particular for the occurrence of druginduced drug resistance. This allows to theoretically justify the clinical notions of drug holiday, i.e., stopping delivering a drug for a defined time interval, and re-challenge, i.e., giving the same drug -and not another one by a wrong consideration that the first one has definitely failed or will fail -after a drug holiday, in cancer [START_REF] Tonini | Rechallenge therapy and treatment holiday: different strategies in management of metastatic colorectal cancer[END_REF].

II. WHERE BIOLOGY MEETS PHILOSOPHY A. Heterogeneity and plasticity of cancer cell populations

Heterogeneity of cancer cell populations [START_REF] Marusyk | Intra-tumour heterogeneity: a looking glass for cancer?[END_REF] has often been mentioned as a major obstacle to the efficacy of drug treatments, in particular of so-called 'targeted therapies' [START_REF] Gillies | Evolutionary dynamics of carcinogenesis and why targeted therapy does not work[END_REF], as treatments that may be very effective on a dominant genetic clone of a tumour cease to be so when another clone emerges [START_REF] Ding | Clonal evolution in relapsed acute myeloid leukaemia revealed by whole genome sequencing[END_REF], be it in malignant haemopathies [START_REF] Ding | Clonal evolution in relapsed acute myeloid leukaemia revealed by whole genome sequencing[END_REF] (aka liquid tumours) or in solid tumours [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF]. This certainly concerns genetic heterogeneity, i.e., situations in which spontaneous mutations have appeared due to the genetic instability of tumour cells, but isogenic heterogeneity resulting from nongenetic instabilities in gene expression due to the lability of epigenetic enzyme activity [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF] is as much important [START_REF] Easwaran | Cancer Epigenetics: Tumor Heterogeneity, Plasticity of Stem-like States, and Drug Resistance[END_REF], [START_REF] Feinberg | Epigenetic modulators, modifiers and mediators in cancer aetiology and progression[END_REF], [START_REF] Mccullough | Plasticity of the neoplastic phenotype in vivo is regulated by epigenetic factors[END_REF], and it occurs before mutations, that remain rare events. The plasticity of cancer cells with respect to the expression of genes [START_REF] Clairambault | Plasticity in cancer cell populations: biology, mathematics and philosophy of cancer[END_REF], [START_REF] Easwaran | Cancer Epigenetics: Tumor Heterogeneity, Plasticity of Stem-like States, and Drug Resistance[END_REF], [START_REF] Mccullough | Plasticity of the neoplastic phenotype in vivo is regulated by epigenetic factors[END_REF], [START_REF] Yuan | Cellular plasticity in cancer[END_REF], as all differentiations, that make an isogenic multicellular organism so varied in cell types (20 for sponges [START_REF] Müller | Review: How was metazoan threshold crossed? The hypothetical Urmetazoa[END_REF], [START_REF] Müller | Origin of the Metazoan Immune System: Identification of the Molecules and their Functions in Sponges[END_REF], about 200 for human beings) is due to epigenetic factors. These are strictly controlled at the organismic level in physiology, and out of control in cancer [START_REF] Easwaran | Cancer Epigenetics: Tumor Heterogeneity, Plasticity of Stem-like States, and Drug Resistance[END_REF], [START_REF] Feinberg | Epigenetic modulators, modifiers and mediators in cancer aetiology and progression[END_REF]. Such plasticity is characterised by the velocity of adaptation of cell phenotypes, such as fecundity, viability or motility, to changes in the cell, tissue and whole organism ecosystem.

B. Philosophy of cancer

A recently autonomised field of knowledge, the so-called philosophy of cancer [START_REF] Bertolaso | Philosophy of cancer. A dynamic and relational view[END_REF], [START_REF] Clairambault | Plasticity in cancer cell populations: biology, mathematics and philosophy of cancer[END_REF], [START_REF] Strauss | Rethinking cancer: a new paradigm for the postgenomics era[END_REF] considers cancer as an impairment of organismic control on both proliferation and differentiations [START_REF] Bertolaso | Philosophy of cancer. A dynamic and relational view[END_REF]. It should be stressed that cancer is a disease that can be met only in the animal kingdom, constituted of multicellular organisms, that, different from plants (plants may also develop tumours, but they remain localised in capsules, never threatening the existence of the organism), are endowed with other properties such as motility and heterotrophy. Multicellular animals, also called Metazoa, emerge from a unique cell, the zygote (fecundated egg) and are thus completely isogenic. The same genome thus exists in all cells of a given animal organism, and the extreme diversity of cell types that constitute an animal is only due to successive differentiations, by selective expression or repression of genes at each differentiation step. Such expressions or repressions are normally under strict epigenetic control, i.e., under the control of epigenetic enzymes that may become impaired in cancer. Furthermore, in cancer, these epigenetic mechanisms are easily manipulated in the plastic cancer cell populations in response to changes of the tumour microenvironment (such as due to the introduction of life-threatening drugs) to develop (reversible) resistant cell phenotypes. The question of understanding what makes a multicellular individual coherent, contrary to a tumour, may not always be summed up by its zygotic origin, and it is the object of the philosophy of immunology [START_REF] Pradeu | The limits of the Self[END_REF], [START_REF] Pradeu | Philosophy of immunology[END_REF], as one of the components that define a multicellular individual is the necessary existence of an immune system [START_REF] Müller | Review: How was metazoan threshold crossed? The hypothetical Urmetazoa[END_REF], [START_REF] Müller | Origin of the Metazoan Immune System: Identification of the Molecules and their Functions in Sponges[END_REF], absent in tumours, which makes possible the success of immunotherapies.

C. The atavistic theory of cancer

In a completely independent way from the above mentioned works, the atavistic theory of cancer was speculated and popularised from 2011 on in a series of papers [START_REF] Davies | Cancer tumors as Metazoa 1.0: tapping genes of ancient ancestors[END_REF], [START_REF] Vincent | Cancer: a de-repression of a default survival program common to all cells?: a life-history perspective on the nature of cancer[END_REF], [START_REF] Lineweaver | Targeting cancer's weaknesses (not its strengths): Therapeutic strategies suggested by the atavistic model[END_REF], relying or partly confirmed in previous [START_REF] Israel | Tumour progression: random mutations or an integrated survival response to cellular stress conserved from unicellular organisms?[END_REF], or following papers examining phylostratigraphic trees [START_REF] Domazet-Lošo | An ancient evolutionary origin of genes associated with human genetic diseases[END_REF], [START_REF] Domazet-Lošo | Phylostratigraphic tracking of cancer genes suggests a link to the emergence of multicellularity in Metazoa[END_REF], [START_REF] Lineweaver | Comparison of the atavistic model of cancer to somatic mutation theory: Phylostratigraphic analyses support the atavistic model[END_REF], [START_REF] Lineweaver | Cancer progression as a sequence of atavistic reversions[END_REF] or archeogenetic data [START_REF] Bussey | Ancestral gene regulatory networks drive cancer[END_REF], [START_REF] Cisneros | Ancient genes establish stress-induced mutation as a hallmark of cancer[END_REF], [START_REF] Trigos | Altered interactions between unicellular and multicellular genes drive hallmarks of transformation in a diverse range of solid tumors[END_REF]. It proposes that cancer is a reversal in evolution from a coherent multicellular organism towards an incoherent, nonviable state closer to unicellularity, but representing an autonomous, coarse and unsuccessful try towards multicellularity from partly organised collections of cells that escape organismic control, a state which these authors name Metazoa 1.0. The interest of considering this atavistic theory is that all cancers are constituted of cells in a same multicellular organism, that have kept in the memory of their common genome defence mechanisms against external aggressions from the environment, that can be unleashed and recruited due to the absence of organismic control, that otherwise would thwart them as introducing too loose cohesion of the ensemble, challenging the very existence of the organism [START_REF] Aktipis | Cancer across the tree of life: cooperation and cheating in multicellularity[END_REF]. Note however that such mechanisms, based on cell de-differentiations and re-differentiations, are at work, strictly controlled, in early development and in adult physiological wound healing.

III. TUMOUR GROWTH AND ITS CONTROL

A. Pitfalls in clinical oncology

In the clinic of cancers, drug treatments meet two main pitfalls: emergence of resistance in cancer cell populations and existence of adverse toxic effects in healthy cell populations. Although most mathematical optimisation works consider one or the other pitfall, it will be shown in the sequel that it is possible to manage them both simultaneously, and not only by following absolute limitations to both population densities, supposed to be taken from clinical considerations, but rather by following with time their respective parts and permanently controlling them.

B. Continuous models of tumour growth

Most models of tumour growth presented until recently were based on ODEs, or PDEs structured in Cartesian space, not many of them are PDEs structured according to phenotypes [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], even fewer according to both phenotypes and space variables [START_REF] Lorz | Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors[END_REF]. The privileged focus of this paper is on structuring tumour growth models in phenotypes [START_REF] Clairambault | A survey of adaptive cell population dynamics models of emergence of drug resistance in cancer, and open questions about evolution and cancer[END_REF].

C. Theoretical therapeutic control

Many models of control, including optimal control methods, mimicking therapeutics in the clinic have been proposed, for ODE or PDE models, and useful references on this topic are in particular [START_REF] Jarrett | Optimal Control Theory for Personalized Therapeutic Regimens in Oncology: Background, History, Challenges, and Opportunities[END_REF], [START_REF] Ledzewicz | On the role of tumor heterogeneity for optimal cancer chemotherapy[END_REF], [START_REF] Schättler | Optimal Control for Mathematical Models of Cancer Therapies[END_REF]. A particular aspect is control in dormancy [START_REF] Djema | Control in Dormancy or Eradication of Cancer Stem Cells: Mathematical Modeling and Stability Issues[END_REF]. In the present article is presented an example of optimal control for a phenotype-structured integro-differential model of the nonlocal Lotka-Volterra type [START_REF] Almeida | Evolution of cancer cell populations under cytotoxic therapy and treatment optimisation: insight from a phenotypestructured model[END_REF], [START_REF] Kuznetsov | Improving cancer treatments via dynamical biophysical models[END_REF], [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF].

IV. MODELS OF ADAPTIVE DYNAMICS

A. Structured population models and adaptive dynamics

The general idea of structured models is to describe heterogeneity in a population of individuals represented by continuous densities by using continuous traits that are relevant to the question at stake, and that are susceptible to evolve with time by adaptation of such traits to changing environmental conditions or interactions with other populations of individuals. This way of representation originally comes mostly from the seminal paper of McKendrick [START_REF] Mckendrick | Applications of mathematics to medical problems[END_REF] on public health problems (see in the original paper the equation p. 122) and from theoretical ecology [START_REF] Diekmann | A beginner's guide to adaptive dynamics[END_REF], however it has been adapted to cell populations, see, e.g., [START_REF] Perthame | Transport equations in biology[END_REF], and a number of following papers. The structuring trait in the population of individuals may be age, size, or physiological traits determining adaptation to the environment such as fecundity, viability [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF], [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF], [START_REF] Lorenzi | Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations[END_REF], motility or plasticity, as will be shown in the sequel. Indeed, as for general and long-term evolution, as regards mutations [START_REF] Michod | Life-history evolution and the origin of multicellularity[END_REF], [START_REF] Nedelcu | Stress responses co-opted for specialized cell types during the early evolution of multicellularity[END_REF], [START_REF] Wagner | Stress-induced evolutionary innovation: a mechanism for the origin of cell types[END_REF], adaptation to stress induced by the environment, as regards epimutations for short-term non-genetic fast phenotypic changes, should be considered as a major determinant of trait differentiation.

B. Tumours and their micro-environment

Populations of cancer cells, in malignant haemopathies or in solid tumours, are amenable to such representations. Their micro-environment may be constituted of other cell populations, in mutualistic [START_REF] Nguyen | Adaptive dynamics of hematopoietic stem cells and their supporting stroma: A model and mathematical analysis[END_REF] or competitive (in particular tumour-immune, for which many models have been published) interactions. It may also consist of mechanical or chemical constraints, as in particular by the introduction of anticancer drugs, that target different aspects of tumour growth, in particular cell death, intrinsic proliferation or action on the structuring traits. A trait of expression of resistance to a given drug has often been considered [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF], [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF], [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], [START_REF] Lorz | Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors[END_REF], [START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF], [START_REF] Pisco | Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: 'What does not kill me strengthens me[END_REF].

C. Tumour strategies to dodge therapeutic control

As mentioned in the introduction, at least two of these strategies have been identified. In the first one [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF], a phenotype of expression of drug resistance structures two cell populations, one tumoral, the other healthy (here healthy only meaning targets of unwanted toxic side effects of the drug), that interact only through a nonlocal logistic term, i.e., an environment integral term I(t), weighted sum of the total mass of the two cell populations, that impinges on the natural death coefficient d(x), representing an assumed competition between all cells for space and nutrients. The model illustrates drug-induced drug resistance, present in cancer, not in healthy cell populations. It will be developed in the next subsection.

The other tumour strategy, that is explored in [START_REF] Alvarez | Evolution of a structured cell population endowed with plasticity of traits under constraints on and between the traits[END_REF] (without control), is related to tumour bet hedging [START_REF] Brutovsky | Structure of Intratumor Heterogeneity: Is Cancer Hedging Its Bets?[END_REF], [START_REF] Gravenmier | Adaptation to Stochastic Temporal Variations in Intratumoral Blood Flow: The Warburg Effect as a Bet Hedging Strategy[END_REF], [START_REF] Philippi | Hedging one's evolutionary bets, revisited[END_REF]. It hypothesises that in order to adapt to unpredictable changes in the tumour environment, tumour cell populations are amenable to display different cell phenotypes, with differentiated abilities to resist different threats. This is represented in the simplest possible way by phenotypic divergence, i.e., the division of the cell population into two phenotypically distinct subpopulations, each one of them maximising a form of fitness related to two different phenotypes related by a flexible (i.e., not simply y = 1 -x) trade-off condition between them.

The main equation, that describes the evolution of a cell population of density n(t, z) at time t bearing phenotype z, where z = (x, y, θ), with for instance x=viability, y=fecundity, and θ=plasticity, under some environmental pressure represented by an advection term V (t, z), runs as:

∂ t n + ∇ • {V n -A(θ)∇n} = (r(z) -d(z)ρ(t))n,
where (V n -A(θ)∇n) • n = 0 for all z ∈ ∂D (n being the outer normal vector on ∂D) and n(0, z) = n 0 (z) for all z ∈ D = {C(x, y) ≤ K} × [0, 1], D defining a trade-off between traits x and y in [0, 1] (e.g., (x -1) 2 + (y -1) 2 ≥ 1), θ free in [0, 1], and the diffusion matrix

A(θ) =   a 11 (θ) 0 0 0 a 22 (θ) 0 0 0 a 33   ,
(where a 11 and a 22 are non decreasing functions of θ) is a representation of how the internal plasticity trait θ impacts the non-genetic instability (due to epigenetic epimutations) of traits x and y, by tuning the diffusion term ∇.{A(θ)∇n}; the advection term

∇.{V (t, z)n} = ∇.{(V 1 (t, z), V 2 (t, z), V 3 (t, z))n},
where for instance V (t, z) = 10 -3 θ(-y, -x, -(x + y)) in simulations, represents the impact of external evolutionary pressure on the population, by changes in the ecosystem hosting the cell population; the integral ρ(t) = D n(t, z)dz stands for the total mass of cells in the population at time t.

The reader is sent to [START_REF] Alvarez | Evolution of a structured cell population endowed with plasticity of traits under constraints on and between the traits[END_REF] to check that phenotypic divergence around two different concentrated phenotypes (x 1 , y 1 ) and (x 2 , y 2 ) is shown to actually occur in the model, while the plasticity phenotype θ simultaneously decreases. In this setting where phenotype space is 3-dimensional, the plasticity trait θ tunes the Laplacian term, while the advection function represents the influence of the environment on the adaptable structure variable (3-d phenotype) z = (x, y, θ) present in the cell population with density n(t, z).

D. A phenotype-structured model with optimal control

The model presented in [START_REF] Almeida | Evolution of cancer cell populations under cytotoxic therapy and treatment optimisation: insight from a phenotypestructured model[END_REF] is a consistent illustration of a drug resistance phenotype-structured model in a cancer cell population with optimal control for the optimised delivery of an anticancer chemotherapy acting as an added death term in a reaction-diffusion setting. However, the somewhat richer integro-differential model presented in [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF] takes simultaneously into account both emergence of drug-induced drug resistance (to be limited) and adverse toxic effects on a nonresistant healthy cell population. As mentioned in the previous section, the two populations are here linked by mixed population terms I H (t) and I C (t), standing for the influence of a global cell environment in the nonlocal logistic term of the Lotka-Volterra model which is described below.

The phenotype x ∈ [0, 1] of expression of drug resistance (x = 0: total drug sensitivity; x = 1: total drug resistance) structures the two populations. However, the functions r H , d H , µ H for proliferation, death and drug sensitivity, respectively, in healthy cells, are chosen much less sensitive to the drug levels and to the environment variables than their correspondent functions r C , d C , µ C in cancer cells.

(Healthy cell population H) ∂ ∂t n H (t, x) = r H (x) 1 + k H u 2 (t) -d H (x)I H (t) -u 1 (t)µ H (x) n H (t, x).
(Cancer cell population C)

∂ ∂t n C (t, x) = r C (x) 1 + k C u 2 (t) -d C (x)I C (t) -u 1 (t)µ C (x) n C (t, x).
Environment terms involving the two cell populations: Each of these equations for the two cell populations is a transport equation of the nonlocal Lotka-Volterra type, structured by the continuous phenotype x of drug resistance, in which the proliferation term r(x) is tuned by a control term representing the effect of a cytostatic (non-cell killing) drug u 2 , while an added death term µ(x) is tuned by the effect of a cytotoxic (cell killing) drug u 1 . As mentioned earlier, the drug-free death term d(x)I(t) stands for global between-cell competition for space and nutrients. The analysis of the asymptotic behaviour of the system in the absence of drugs (or with constant drug doses) shows concentration of the phenotype x in both populations on a discrete set in [0, 1] [49], which in simulations appears as reduced to a unique point. Then comes the question of optimally controlling this asymptotic behaviour by a combination of the two time-continuous (or L ∞ ) drugs (u 1 , u 2 ).

I H (t) = a HH .ρ H (t) + a HC .ρ C (t), I C (t) = a CH .ρ H (t) + a CC .ρ C (t), with ρ H (t) = 1 0 n H (t, x) dx, ρ C (t) =
Optimal control problem OCP: find controls (u 1 , u 2 ) minimising in fixed horizon T the objective function

C T (u 1 , u 2 ) = ρ C (T ) = 1 0 n C (T, x) dx under the constraints ρ H (t) ρ H (t) + ρ C (t) ≥ θ HC , ρ H (t) ≥ θ H .ρ H (0)
(the last constraint, with, e.g., θ H = 0.6, to limit damage to healthy cells) Then it can be shown that the optimal therapeutic trajectory (u 1 , u 2 ) in large time T > 0 consists of 3 parts:

• a long-time part, with constant controls on [0, T 1 ] (where

u 1 = 0, u 2 saturated along the constraint ρ H (t) ρ H (t) + ρ C (t)
= θ HC ), at the end of which populations have concentrated in phenotype (for T 1 large); • a free arc (no constraint saturating) where controls

u 1 = u max 1 and u 2 = u max 2 ; • a boundary arc along the constraint ρ H (t) ≥ θ H .ρ H (0) with u 1 constant (< u max 1
) and u 2 = u max 2 . The reader is sent to [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF] for proofs, numerical simulations and illustrations. The figure below shows the resulting optimal control (u 1 , u 2 ) for θ HC = 0.4 and θ H = 0.6.

V. DISCUSSION ABOUT STRUCTURED MODELS A. Including targets for control

In all mathematical models, be they phenotype-structured or not, that are designed in the perspective of therapies, Fig. 1. Evolution with time t of the x resistance phenotype distribution curves n H (t, x) and n C (t, x) and solution (u 1 , u 2 ) of the OCP problem for T = 30. In the lower two panels, one can see that the cytotoxic drug u 1 level is maintained positive, but as low as possible, close to zero, for a long time, before delivering its maximal tolerated dose on a short time interval, so as to keep the cell population as sensitive as possible until time T 1 (here T 1 ≈ 20) -which is the drug holiday strategy used in the clinic of cancers [START_REF] Tonini | Rechallenge therapy and treatment holiday: different strategies in management of metastatic colorectal cancer[END_REF] -, then delivered at a medium dose. The cytostatic drug u 2 is maintained at a low level, but not close to zero, until some time before and close to T 1 , and then delivered at its maximum tolerated dose. One can see on the two rightmost panels that the cancer cell population ρ C goes to extinction while the healthy cell population ρ H , target of unwanted toxic side effects of the drugs, is preserved above a prescribed level. The central panel shows the evolution of the ratio of drugsensitive cancer cells ρ CS (t) = 1 0 (1 -x)n C (t, x)dx over total cancer cells ρ C (t) = 1 0 n C (t, x)dx. Performed in the modelling language AMPL with expert simulation routines IpOpt. Adapted from [START_REF] Pouchol | Asymptotic analysis and optimal control of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy[END_REF] with permission.

physiologically based targets of control must be included in the equations. If the control primarily targets cell death (case of most classic chemotherapies), then it must be present in the equations as an added death term. If it targets proliferation, as for, e.g., tyrosine kinase inhibitors (TKIs) that are also antagonists of epidermal growth factor receptors (EGFR antagonists), then it must be present as tuning the intrinsic (that is, defined by setting death terms to zero) instantaneous growth term ∂n n∂t . If it is an epigenetic drug that controls nongenetic instability, then it must address a term tuning a Laplacian ∂ 2 n ∂x 2 or other terms representing nongenetic instability in a reaction-diffusion equation. Controlling the tumour microenvironment might also be thought of as tuning an added advection term ∂n ∂x , the structure variable x being here assumed to be one-dimensional for simplicity of presentation.

B. Taking into account heterogeneity

Various ODE models claim to take account of cell population heterogeneity by compartments representing discrete heterogeneity, such as (totally) drug-sensitive and (totally) drug-resistant cell subpopulations. This relies on a binary vision of structuring phenotypes, exchanges between compartments representing switching between discrete subpopulations, sensitive vs. resistant, as in the SIR models for infectious diseases introduced by Kermack and McKendrick in 1927 [START_REF] Kermack | Contribution to the mathematical theory of epidemics[END_REF]. The vision presented in this paper is closer to the type of transport equation introduced by McKendrick (alone) in 1926 [START_REF] Mckendrick | Applications of mathematics to medical problems[END_REF], to which nonlocal terms may be added, with or without diffusion and advection terms. Indeed, there is no reason why functional structuring traits in the population, such as fecundity, viability, motility and plasticity should be discrete. This is the reason why heterogeneity should a priori be represented by continuous traits. It depends on the underlying questions at stake, and on the ways of appreciating such heterogeneity. If the available measures and if the biological questions may be summed up by discrete characters, why not design compartmental ODE models to take account of heterogeneity? However, in general, models structured in continuous traits, yielding integro-differential equations, or partial differential equations with diffusion and advection terms, should be the rule, as heterogeneity in biology is only coarsely represented by discrete variables.

C. Plasticity and tumour microenvironment

As mentioned earlier, plasticity understood as nongenetic instability may be represented by a Laplacian. This is the case in [START_REF] Alvarez | Evolution of a structured cell population endowed with plasticity of traits under constraints on and between the traits[END_REF], [START_REF] Chisholm | Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation[END_REF], [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation[END_REF], where nongenetic instability or plasticity play a major role in the evolution of the solutions. The influence of evolutionary pressure (or stress) due to the tumour micro-environment is adequately represented in both cases by an advection (or drift) term, the specific effect of which is studied in [START_REF] Chisholm | Effects of an advection term in nonlocal Lotka-Volterra equations[END_REF], and interpreted as displacing the expected concentration point of the one-dimensional structure variable of the nonlocal Lotka-Volterra reactiondiffusion under consideration. In [START_REF] Alvarez | Evolution of a structured cell population endowed with plasticity of traits under constraints on and between the traits[END_REF], plasticity is endowed with the nature of an intrinsic trait tuning the Laplacian and evolving with time as the other (competing) traits. It is noteworthy that in order to represent competing phenotypes, such as fecundity and viability, a dimensionality of the space of phenotypes should be greater than one. If we add Cartesian space to phenotype space (as performed in [START_REF] Lorz | Modeling the effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors[END_REF], however with only one dimension in each case), this will make the model even more complex to study. Nevertheless, this is the framework of development of actual solid tumours, in which local (in space) fixation of phenotypes may lead with genetic instability to the emergence of spatially isolated mutations according to a Darwinian-like phylogenetic tree as has been shown in [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF]. In other words, mixed space and phenotype models should consider firstly concentration of phenotypes, reversible and costly in terms of energy, likely using epigenetic enzymes to face micro-environmental stress [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], in spatially isolated regions (spatial basins of attraction for such concentration), and secondly, under the pressure of microenvironmental stress in these regions, possible emergence of mutations fixing -on a less costly mode than through the action of epigenetic enzymes -a favourable phenotype.

D. Interacting structured cell populations

Interactions of tumour cell populations with surrounding healthy ones may be mutualistic [START_REF] Nguyen | Adaptive dynamics of hematopoietic stem cells and their supporting stroma: A model and mathematical analysis[END_REF] or competitive. Structured models of tumour immune-interactions, for the latter, are in development [START_REF] Kaid | A phenotype-structured model for the tumour-immune response[END_REF], since the beginning of modern immunotherapies, and they offer many opportunities, including for the development of control methods, to mathematicians.

VI. CONCLUSION AND PERSPECTIVES

In this short review of some mathematical models of tumour growth and its control, focused on structured cell population models, principles of physiologically based modelling for cancer growth, to take into account cell population heterogeneity and plasticity, with examples of models with and without control, have been presented. Clearly, the categories of "cancer" and "healthy" cell populations we have considered in our model with optimal control are very general. The illustrated solution to the above mentioned optimal control problem OCP proposes a qualitative strategy of control that should be adapted to defined clinical settings, in order to be actually useful in practical oncology. Nevertheless, the proposal is on the table, awaiting for therapists to consider it. In the immediate future, modelling tumour-immune interactions with representation of modern immunotherapies involving tumour control by immune checkpoint inhibitors (ICIs, that boost immune cells) is another goal which is presently pursued [START_REF] Kaid | A phenotype-structured model for the tumour-immune response[END_REF], with the same principle of structuring cell populations according to a continuous phenotype representing their heterogeneity. However, even immunotherapies act by tumour cell kill induced by drugs, and their adverse events, although rare, may be severe, even fatal, as in the case of cardiac or diaphragmatic myopathies. The grail of anticancer treatments might be, in particular by using redifferentiation therapies, to avoid cancer cell kill and its adverse toxic events on healthy tissues. Unfortunately, this has been achieved so far only in very special cases, such as acute promyelocytic leukaemia (APL) or chronic myelogenous leukaemia (CML) with drugs All-TransRetinoic Acid (ATRA) and Imatinib Mesylate, respectively. Progression in the understanding of mechanisms of tissue cohesion in multicellular organisms and their recovery after local impairment in cancer might offer new opportunities for therapeutics, although most of these mechanisms are still widely unknown.

1 0n

 1 C (t, x) dx, and u 1 cytotoxic, u 2 cytostatic drug concentrations.


Valentin Boussot
email: valentin.boussot@univ-rennes1.fr

Jean-Louis Dillenseger
email: jean-louis.dillenseger@univ-rennes1.fr

Statistical model for the prediction of lung deformation during video-assisted thoracoscopic surgery

Keywords: General statistical motion model, Images Registration, CBCT, Lung deformation, Video-assisted thoracic surgery (VATS)

Accurate lung nodule localization during video-assisted thoracic surgery (VATS) for the treatment of early stage lung cancer is a surgical challenge. Recently, a new minimally invasive approach for nodule localization during VATS has been proposed, which consists in compensating by a biomechanical model the very large lung deformations occurring before and during surgery. This estimation of the deformations allows to transfer the position of the nodule visible on the preoperative CT to an acquisition of the lung performed during the operation using a Cone-Beam CT scanner (CBCT 2). But, in this approach, an additional CBCT acquisition (CBCT 1) must also be acquired just after the patient is placed in the operative position in order to estimate the deformations due to the change of the patient's position, from supine during the CT acquisition to lateral decubitus in the operating room. Our goal is to simplify this procedure and thus reduce the radiation dose to the patient. To this end, we propose to improve this solution by replacing the CBCT 1 acquisition by a model allowing to predict these deformations. This model is defined using the lung state information from CBCT 2 and a general statistical motion model built from the position change deformations already observed on other patients. We have data from 17 patients. The method is evaluated with a leave-one-out cross-validation on its ability to reproduce the observed deformations. The method reduces the average prediction error from 12.14 mm without prediction to 8.12 mm for an average prediction, and finally to 6.33 mm for a prediction with our model fitted to CBCT 2 only.

INTRODUCTION

Lung cancer is the leading cause of cancer death in the world for both men and women. Its mortality rate depends largely on the stage at which it is diagnosed. In the case of early detection, i.e. when cancer cells are still confined to the tissue of origin, the reference treatment is video-assisted thoracic surgery (VATS). The operation consists of removing the small area of the lung that contains the lesion under video guidance in a procedure called Wedge. To perform this procedure, the surgeon makes 3 incisions of about 3 cm in diameter in the chest, one for the endoscope and two to introduce the instruments. This minimally invasive operation allows a better and faster recovery and gives the patient a good chance of survival.

The localization of the pulmonary nodule during this operation is a prerequisite for its successful completion. However, the precise localization of the nodule to within 5 mm presents a difficulty from a surgical point of view. Indeed, the nodule is not visible from the surface of the lung and the small incisions to access the lung do not allow the surgeon to palpate it. The nodule is however visible on the preoperative CT, but there is a great anatomical disparity between the preoperative and intraoperative configurations, due to the large pulmonary deformations introduced during surgery. These deformations are of two kinds. On the one hand, the patient is in dorsal decubitus position during the CT, but he is placed in lateral decubitus on the operating table. On the other hand, the incisions made by the surgeon allow air to enter the pleuropulmonary space which produces a strong collapse of the lung into the mediastinum called pneumothorax. These deformations greatly modify the shape and internal pulmonary structures and thus the position of the nodule in relation to the reference CT.

In current practice, localization is based on the introduction of a marking device (needles, hook wires, microcoils or dyes) into the lung during an initial preoperative percutaneous procedure under fluoroscopic and CT guidance. However, this invasive percutaneous preoperative procedure leads to high complication rates. In order to avoid this preoperative procedure, Rouzé et al proposed a totally non-invasive localization approach based on intraoperative imaging only. [START_REF] Rouzé | Small pulmonary nodule localization with cone beam computed tomography during video-assisted thoracic surgery: a feasibility study[END_REF] A Cone-Beam CT scan (CBCT 2) is performed during the operation. This acquisition allows first to identify and localize the nodule in the CBCT, and then to guide the surgeon to this nodule. However, this solution suffers from a problem of reliability in the case of low density or deep nodules. Indeed, on the one hand, the quality of the CBCT image is reduced compared to a CT image with a lower contrast between the nodule and the surrounding lung parenchyma and on the other hand, pneumothorax leads to atelectasis, i.e. a high densification of the lung tissue. These two problems make the segmentation of these types of nodule impossible. To overcome this problem, Alvarez et al proposed to model the deformations of the lung by a biomechanical model to transfer the position of the nodule visible in CT on the CBCT 2. [START_REF] Alvarez | A hybrid, imagebased and biomechanics-based registration approach to markerless intraoperative nodule localization during video-assisted thoracoscopic surgery[END_REF] As the deformations of the lung are very large and of different origins, Alvarez et al in their approach decided to acquire an intermediate CBCT volume (CBCT 1) just after the patient is placed in the operative position but before the pneumothorax. This CBCT 1 allows: 1) to estimate the deformations due to the change in the patient's position by a registration between the preoperative CT and the CBCT 1 and 2) to establish the initial position of the biomechanical model in order to simulate the deflation (pneumothorax) of the lung towards the surfaces seen on the CBCT 2. This estimation is essential for a good simulation of the pneumothorax and its allows to compensate initially for non-rigid deformations of the order of 1 cm on average and up to 30 mm for some clinical cases.

The objective of our work is to simplify this protocol and to reduce the radiation dose per patient by eliminating the intermediate CBCT 1 acquisition. To do so, instead of this acquisition, we will try to predict the deformations due to the change of the patient's position from the CT and CBCT 2 alone. This prediction will be based on a statistical model of movements established from the deformations already observed on the data of the previous patients who had the 3 acquisitions (CT, CBCT 1 and CBCT 2).

MATERIALS AND METHODS

Data

We have data from 17 patients, 12 patients included for a wedge nodule resection performed on the right lung and 5 patients for a nodule present in the left lung. The VATS surgical intervention was performed at the Rennes University Hospital, Rennes, France. This study was approved by the local ethics committee and the patients gave their informed consent prior to the procedure.

For each patient we have three volumes: 1) a preoperative CT, where the patient is in the supine position and was instructed to hold his breath during the capture (end of inspiration cycle); 2) a CBCT acquisition, called CBCT 1, performed during the operative phase, i.e. the patient is in lateral decubitus under sedation and mechanical ventilation. However, this acquisition is performed before the incisions therefore the lung is completely inflated; and 3) a CBCT acquisition, called CBCT 2, made after the thorax perforation with the resulting pneumothorax. Figure 1 shows an axial view of these 3 images.

We have built one deformation model for the right lung and one for the left lung. However, for the sake of clarity, in this paper we will present and evaluate only the deformation model for the right lung. The construction of the model for the left lung followed the same methodological approach.

Method workflow

The global workflow of the study is presented Figure 2.

In the learning phase, we had at our disposal the CT-CBCT 1 image pair for each patient. For our protocol we need the lung lobes segmented in CT and the spine segmented in CBCT modality. The idea is: a) for a given patient, to estimate the deformation field of the change of pose between CT and CBCT 1. For this purpose, CBCT 1 is aligned to the CT by a rigid spine registration. Then, the deformation T (.) of the lung due to the change of pose is estimated by elastic registration between CT and CBCT 1; b) to establish a matching field between the lungs of the several patients. For this, we align the CT of the different patients to a model lung. This is done by inter-patient registration of the segmented lungs in the CT. This matching field T i,t (.) between each patient's lung and the model lung ensures that each sampled vector accounts for the deformation of the same lung region regardless of the patient; c) create a model of lung tissue deformations using a weighting PCA computed on a parsimonious sample x j of the different deformation fields describing the pose change. This allows us to create a deformation model with degrees of freedom corresponding to the main deformation modes observed.

In the prediction phase for a new patient, we have at our disposal the preoperative CT and the CBCT 2 after pneumothorax. Our starting assumption is that the pleuropulmonary space does not undergo deformation between CBCT 1 and CBCT 2 but undergoes significant deformation by the change in patient position. The segmented pleuropulmonary space on CBCT 2 thus gives the shape of the lung surface after deformations due to the change of pose. The prediction of the change of pose deformation can be formulated as an inverse problem. For this purpose: a) CBCT 2 is aligned to CT by rigid spine registration; b) CT is aligned to the model lung by inter-patient registration; c) the principal deformation modes are identified by an inverse problem from the lung surface in CT and the pleuropulmonary space surface in CBCT 2; d) once the principal modes are identified, the model vectors are associated with a deformation field by a third order B-Spline interpolation.

The key points of our method (segmentation, registration, creation of the statistical model of deformations and resolution of the inverse problem) will now be detailed.

Segmentations

In our method we need to segment 3 anatomical structures: 1) the 3 lobes of the right lung or the 2 lobes of the left lung in the CT modality, 2) the spine on the CBCT 1 and CBCT 2 images and 3) the pleuropulmonary space visible in the CBCT 2.

1) The anatomical description of the lobes was performed with the open source tool Chest Imaging Platform (CIP) present on the 3D Slicer software. [START_REF] Jose Estepar | Chest imaging platform: an open-source library and workstation for quantitative chest imaging[END_REF] It is a semi-automatic tool that segments the 5 lobes of both lungs. This algorithm must be manually initialized with a set of markers to help identify the lobe boundaries.

2) In our case, a coarse segmentation of the spine is sufficient, it is performed by an automatic algorithm based on the intensities with a priori on the position of the spine that depends on the lung that is operated.

3) The description of the pleuropulmonary space in CBCT 2 is the union of the deflated lung and the air present in the pleural space. This description was achieved by an iterative region growing approach implemented by the 3D slicer software. The implemented algorithm requires manual initialization and some interactive postprocessing.

Large deformation image registration method

In our method, two steps require the estimation of a large deformation field between two images: the estimation of the deformation fields between CT-CBCT 1 for one patient and the inter-patient alignment of the CTs. To solve them we used the classical multiresolution iterative registration method based on image intensity implemented in Elastix. [START_REF] Marstal | SimpleElastix: a user-friendly, multi-lingual library for medical image registration[END_REF] In both cases, the similarity measures measured between the fixed image I f and the moving image I m 1) uses a 3rd order B-Splines β 3 (.) deformation model with control points x k arranged on a homogeneous grid with σ = 20 mm spacing on the highest resolution in the fixed image domain Ω f .

T µ (x) = x + x k ∈Ω f p k β 3 ( x -x k σ ) µ = {p 1x , p 1y , p 1z , ..., p N x , p N y , p N y } (1)

CT-CBCT 1 registration

The estimation of the deformations due to the change of patient's pose is performed by an elastic registration of the CT-CBCT 1 image pair. This estimation is done in 2 steps, 1) a rigid registration and 2) an elastic registration.

1) The rigid registration aims at aligning the structures not affected by the deformations in order to estimate in a second time the mean deformation and the main deformation modes with deformation fields relative to this fixed reference point. For this purpose, the two volumes are initially aligned to the spine. This is done by first applying a 90 degree transformation, and then the two volumes are aligned using two matched markers manually placed on one of the vertebrae in the CT and CBCT 1. A rigid registration is then performed using a normalized cross-correlation similarity measure estimated from 3000 samples on the spine in the CBCT.

2 ) One of the key points of the elastic registration is the cost function that estimates the displacement of control points between the fixed image I f and the moving image I m . In our case, we decided to integrate two similarity measures in our cost function: a normalized cross-correlation (N CC) on the amplitudes of the image gradients as originally proposed by Haber et al [START_REF] Haber | Intensity Gradient Based Registration and Fusion of Multi-modal Images[END_REF] and a normalized mutual information on the image intensities (M I):

C(µ; I f , I m ) = -N CC(µ; ||∆I f || 2 , ||∆I m || 2 ) -αM I(µ; I f , I m ), α = 1 (2)
where ||∆I|| 2 denoting a magnitude operator on the gradient of the image I.

The two similarity measures are estimated from 1000 samples extracted within the CT lung volume. The few degrees of freedom of the deformation model allow to have a good description of the deformation by posing the problem well enough to avoid the need of regularization.

The choice of hyperparameters for these algorithms was found by an exhaustive search using a subset of patients (6 clinical cases) for whom markers matched on both images were placed by the pulmonary surgeon (≈ 30 for each case). We chose the combination that yielded the lowest Target Registration Error (TRE) between the markers. The median error was less than 1 mm, which is consistent with the resolution of the voxels in the images

Inter-patient lung matching

The matching field between each patient's lung and a model lung is achieved by the registration of each of the lobes of both lungs. The model lung was selected as the lung with the largest volume in the database. This choice is important because the displacement vectors will be sampled spatially in the model lung and a contraction the spacing of the sampled positions rather than a dilation is desirable.

The lungs between patients show large disparities in size and shape. The challenge of lung registration is not to match their fine structures to each other but to match anatomically close regions. To do this, an affine registration process first compensates for differences in size and orientation. A non-rigid registration is then applied to compensate for the difference in shape. As we want to perform a lobe registration without influence from the outside, as Wu et al [START_REF] Wu | Evaluation of deformable registration of patient lung 4DCT with sub-anatomical region segmentations[END_REF] we assign a very low value (-2000) to the voxels outside the lobes. The similarity measure to be minimized is the combination of the cross-correlation (CC) and a penalty term P (bending energy) that penalizes sharp changes in local transformation curvature in order to make the deformations as smooth as possible while fitting the shape of the lung model. The cost function takes the following form:

C(µ; I f , I m ) = -CC(µ; I f , I m ) + αP (µ) α = 10 and P (µ) = 1 N N i ||H (Tµ(xi)) || 2 F ( 3
)
where H (Tµ(x)) is the Hessian matrix of the transformation T µ (.).

Weighting principal component analysis

A PCA model was used to observe the main deformation modes of the lung population from its mean deformation.

To form our statistics, we used two types of sampling. On the one hand, a volume sampling using a grid of points in the model lung region with a spacing of 16 mm and which includes 866 points. This type of sampling aims to recover the initial deformation field with a B-Spline interpolation. On the other hand, a sampling of the lung surface of the model lung. These sampling points correspond to the nodes of a triangular mesh obtained by a smooth mesh of the model lung segmentation for a total of 7218 points. This type of sampling allows to associate a lung shape to the model.

For each patient i, the deformation field modeling the pose change is sampled using these points, which are previously matched to the patient's CT lung using the model lung matching field. Each sampled vector is flattened and concatenated into a vector ⃗ x i ∈ R (866+7218) * 3 . The set of vectors in the training base thus formed is assembled into a matrix X with I rows (number of patients) and J × 3 columns (Number of vectors × 3). Formally:

X i,j = T (T i,t (x j )) -T i,t (x j ) (4)
One of our major problems is that CBCT field of view (FOV) is limited and does not allow to capture the whole lung. The estimation of its deformation is only valid in this FOV. In order not to bias the mean deformation and the main observed deformation modes, we defined a weight that quantifies for each patient and for each vector the confidence we have in the described estimation. The weights corresponding to the vector are functions of their arrival position x j in the CBCT 1: in the inside of the FOV the weight is fixed at 1 and in the outside the weight decreases exponentially with respect to its distance to the CBCT FOV with a half-life fixed at a distance of p = 10 mm:

w * i,j = e log(0.5) p d f ov i (xj ) w i,j = w * i,j I k=1 w * k,j (5)
where d f ov i (.) corresponds to a distance map from the FOV of the CBCT 1 after rigid registration of the CT/CBCT 1 for the patient i. It is obtained by a unsigned Danielson distance map 7 calculated from a cylindrical mask adjusted to CBCT 1 and resampled in the CT domain. The weight matrix W is normalized so that the columns are summed to 1.

Figure 3 shows the 2 kind of sample points (volume in yellow and surface in red) and the confidence map quantified by weights (color map). The weighted average deformation is calculated as follows:

µ j = I i=1 w i,j x i,j (6)
The construction of the deformation model poses the problem of the low rank approximation of the matrix X. Let be the following optimization problem:

arg min P,B χ 2 = ij w ij [x i,j -µ j - M k=1 p k,j b k,i ] 2 (7)
with P = (p k,j ) 1≤k≤M,1≤j≤J a matrix whose column are the principal components to find, and B = (b k,i ) 1≤k≤M,1≤i≤I is a matrix of coefficients to fit X using P . The mode number M , is assumed to be small and P is subject to an orthonormality constraint, P T P = 1. To solve this optimization problem we used an implementation of an Expectation-Maximization type algorithm proposed by Bailey (EMPCA). 8

Inverse problem formulation

The deformation modes for a new patient are searched so that the surface of the CT lung deformed by the model fits the pleuropulmonary space surface extracted from CBCT 2. For this purpose, the CBCT 2 is first aligned on the CT using a rigid registration of the spine in order to remove the rigid components. The CT is also aligned to the lung model using the registration described in section 2.4.2.

The objective is then to estimate the deformation modes from the measurement of the distances between the mean deformation surface described by the lung model and the pleuropulmonary space surface of CBCT 2. For a quick measurement of these distances, the surface nodes of the mean model are immersed in an square signed Danielson distance map 7 calculated from the pleuropulmonary space mask resampled in the CT volume domain (see Figure 4). The main mode B is estimated by the following optimization problem: arg min

B L i = 1 N j∈Ωj d(µ j + P T j B) (8)
where d(.) corresponds to a distance map from pleuropulmonary mask and P T j ∈ R 3×M a principal component matrix of the j th vector from a subset of vectors localized at the lung surface in the FOV of CBCT 2. The points on the medial face and the diaphragm region are removed from the set.

To make the optimization method more robust to outliers, we applied a non-linear function that allows to give less importance to these points. To solve this optimization problem we applied adagrad:

B t+1 = B t - η g t g T t g t with B 0 = [0, 0, ...], η = 1 and g t = ∂L ∂B = 1 N i P j d ′ (µ j + P T j B) 2 d(µ j + P T j B) (9)
where d ′ (x) = ∂d(x) ∂x is the gradient image of the distance map.

EVALUATION

Evaluation protocol

The experiments were conducted in a leave-one-patient-out procedure on the 12 patients operated on the right lung. We propose a quantitative evaluation of our method. For a test patient, we predict, using our model, CT and CBCT 2, the pose change deformations. We then compare this prediction to the deformations estimated by the elastic registration of the CBCT 1 to the CT used for the learning as described in the workflow (section 2.2) and considered as reference truth. The criterion consists in measuring on the test data, the Target Registration Error (TRE) which is the L2-norm between the volume sampling points x j projected by the estimated deformation T (x j ) (reference truth) and by the predicted deformation T (x j ) within the FOV of CBCT 1:

T RE = ||T (x j ) -T (x j )|| 2 (10)
We compared the performance obtained by 4 prediction models: (Rigid ) no prediction (rigid registration between CT and CBCT 2); (Mean) the prediction of the mean deformation of our model; (Optimal ) the prediction with the estimation of the deformation modes by solving the inverse problem of our method on the pleuropulmonary surface of CBCT 1; and, (Prediction) the prediction with the estimation of the deformation modes by solving the inverse problem of our method on the pleuropulmonary surface visible in the CBCT 2.

The aim of the (Optimal ) test was to verify the accuracy of the model using the same conditions (pleuropulmonary space) as during training. On the other hand, in clinical routine this prediction model cannot be used because we will no more have access to CBCT 1.

Results

Figure 5 and Figure 6 show the error distribution for all the 12 patients according to the 4 predictions. The average median error for each patient for each prediction is (Rigid ) 12.79 mm, (Mean) 8.09 mm, (Optimal ) 3.86 mm and (Prediction) 6.33 mm.

Figure 7 shows the modeling performance of our wPCA model on the test data as a function of the number of M modes used to fit the training data (see (7)). The performance of the 12 models is estimated in the form of the characteristic values of the distribution of errors on the 12 test patients: median value (50%), first and third quartiles (25% and 75%) and limits (5% and 95%).

DISCUSSION AND CONCLUSION

This paper presents a method for predicting lung deformations resulting from the change in the patient pose between preoperative and intraoperative configurations during a VATS procedure. The problem was approached as an inverse problem with the direct problem being a deformation model inferred from pose change deformations estimated on an intermediate imaging applied to a set of patients.

The accuracy of the deformation prediction was evaluated against its estimation using a leave one out cross validation (LOOCV). A simple rigid registration results in a average TRE median of 12.14 mm. This error is reduced to 8.09 mm by fitting an average model to the patient's data and to 6.33 mm by estimating the correct mode by solving the inverse problem.

These results were obtained with relatively few patients. However, our experiments show that we have some generalization in that, according to Fig. 7, we are able to significantly compensate for distortions on the test data with only 4 modes.

Our initial assumption is that the pleuropulmonary space did not undergo deformation between the change of pose (CBCT 1) and the pneumothorax (CBCT 2). In this case, the mode estimation is accurate and the cost function is convex. The direct model thus seems to pose the problem well, by capturing the link between the deformation of the lung surface and its internal deformation (cf. the results obtained by the (Optimal) mode). However, this hypothesis is not verified in reality. Pneumothorax often results in a collapse of the mediastinum and an upward movement of the diaphragm. Secondly, the costal pleura also undergoes a slight deformation during pneumothorax. This, combined with the problem of the limited field of view of CBCT 2, leads to results that are somewhat less accurate than expected (cf. the results obtained by the (P rediction) mode).

One solution that has been explored to solve these problems has been to incorporate into our model the deformations of the pleuropulmonary space associated with pneumothorax, however these deformations are not regular and in some cases may even degrade the performance of the PCA model. The success of a good solution of the inverse problem relies on the choice of the useful region points. Our experiments have shown that among the 4000 points of the lung surface initially present in the field of view of the CBCT 2, about 2500 are useful and lead to an accurate estimation. The problem now is to find the right heuristics to automatically identify the relevant regions. Future work will consist of improving the inverse problem solving scheme. Moreover, we plan to evaluate the prediction of the deformations of the patient's position change on the complete framework, 2 i.e. on its capacity to initialize well the biomechanical model in order to predict the position of the nodule.

Figure 1 .

 1 Figure 1. Slices in approximately the same transversal plane of, from left to right, the preoperative CT (supine position); the intraoperative CBCT 1 (lateral decubitus position) and intraoperative CBCT 2 with the pneumothorax

Figure 2 .

 2 Figure 2. Workflow of our study

Figure 3 .

 3 Figure 3. CT-model volume with the 2 types of sample points (volume in yellow and surface in red) and the confidence map. The confidence function is encoded by a color map, the purple region corresponds to the CBCT 1 field of view

Figure 4 .

 4 Figure 4. Superimposed on a slice of the CBCT 2 volume: in green the contour of the pleuropulmonary space in CBCT 2 with the distance map to this contour represented as a color map; in yellow the contour of the CT lung aligned on the CBCT 2; in blue the boundary of the mean deformation surface described by the lung model, the red markers are the surface nodes of the mesh projected in the distance map.

Figure 5 .

 5 Figure 5. Prediction performance for the first 6 patients

Figure 6 .Figure 7 .

 67 Figure 6. Prediction performance for the last 6 patients

ACKNOWLEDGMENTS

The work presented in this article was supported by the Région Bretagne through its Allocations de Recherche Doctorale (ARED) framework and by the French National Research Agency (ANR) as part of the VATSop project (ANR-20-CE19-0015).