Aurélien Marc

Romain Marlin

Flora Donati

Mélanie Prague

Marion Kerioui

Cécile Hérate

Marie Alexandre

Nathalie Dereuddre-Bosquet

Julie Bertrand

Vanessa Contreras

Sylvie Behillil

Pauline Maisonnasse

Sylvie Van Der Werf

Roger Le Grand

Jérémie Guedj

Impact of variants of concern on SARS-CoV-2 viral dynamics in non-human primates

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

Introduction

The sever acute respiratory coronavirus 2 (SARS-CoV-2) is the causative agent of the Coronavirusinduced disease 2019 (COVID-19) cumulating more than 500 million cases and over 18 million death as measured by excess mortality as the end of 2022 [START_REF]COVID-19) Dashboard | WHO Coronavirus (COVID-19) Dashboard With Vaccination Data[END_REF][START_REF] Wang | Estimating excess mortality due to the COVID-19 pandemic: a systematic analysis of COVID-19-related mortality, 2020-21[END_REF]. Repeatedly, several variants have emerged and although most of them vanished quickly, some of them, called Variants of Concern (VoC), in particular alpha, beta, gamma, delta and omicron have caused dramatic epidemic rebounds [START_REF] Campbell | Increased transmissibility and global spread of SARS-CoV-2 variants of concern as at June 2021[END_REF][START_REF]COVID-19 weekly epidemiological update[END_REF][START_REF]COVID-19 weekly epidemiological update[END_REF]. These variants have acquired specific mutations enhancing their infectious capacities and escaping the immune response, leading to a dramatic loss of efficacy of monoclonal antibodies [START_REF] Choudhary | Emergence of SARS-CoV-2 Resistance with Monoclonal Antibody Therapy[END_REF]. They have also caused a large drop in vaccine efficacy against disease acquisition even though until now vaccine remain largely effective against severe disease [START_REF]Effect of Covid-19 Vaccination on Transmission of Alpha and Delta Variants[END_REF][START_REF]Virological characteristics of SARS-CoV-2 vaccine breakthrough infections in health care workers | NCRC [Internet[END_REF][START_REF] Andrews | Covid-19 Vaccine Effectiveness against the Omicron (B.1.1.529) Variant[END_REF]. While several millions of individuals have been infected by these VoC, we still do not have a precise understanding on the effects of VoC on viral load. Even though some effects on larger levels of viral excretion have been reported [START_REF] Elie | Variantspecific SARS-CoV-2 within-host kinetics[END_REF][START_REF] Blanquart | Characterisation of vaccine breakthrough infections of SARS-CoV-2 Delta and Alpha variants and within-host viral load dynamics in the community, France, June to[END_REF][START_REF] Naveca | COVID-19 in Amazonas, Brazil, was driven by the persistence of endemic lineages and P.1 emergence[END_REF][START_REF] Li | Viral infection and transmission in a large, welltraced outbreak caused by the SARS-CoV-2 Delta variant[END_REF], these studies often lack of robustness, and may be biased by many confounding factors that complicate comparisons, in particular reporting biases, heterogeneity in the incubation period and vaccination coverage. In that context where human clinical data are difficult to interpret, the non-human primate (NHP) experimental model offers a unique opportunity to describe infection with SARS-CoV-2 in detail in a fully controlled environment. Since 2020, our group has conducted many studies to evaluate the effects of antiviral drugs or vaccines in this model [START_REF] Maisonnasse | Hydroxychloroquine use against SARS-CoV-2 infection in non-human primates[END_REF][START_REF] Maisonnasse | COVA1-18 neutralizing antibody protects against SARS-CoV-2 in three preclinical models[END_REF] , and showed its large predictive value [START_REF] Antonio | SARS-CoV-2 viral dynamics in non-human primates[END_REF]. Here, we analysed retrospectively viral load data obtained in 78 animals that were included as control arms of these studies and that were infected with different strains of SARS-CoV-2 (historical, beta, gamma, delta and omicron (BA.1)). In addition, we performed longitudinal measures of viral culture to evaluate a potential effect of VoC on viral infectivity. Using the techniques of mathematical modelling, we characterize the viral kinetics in these animals and we discuss their biological insights.

Results

Variant of concern viral kinetics

Several biomarkers were measured, both genomic RNA and subgenomic RNA were quantified at regular interval over all the study period and infectious titers at 2 times points. All macaques developed a rapid infection with genomic viral load peaking between 2-and 3-day post-infection (dpi) for the historical and beta variant, 3.5 dpi for variant delta and 4 dpi for variants gamma and omicron (BA.1). Genomic viral load was cleared at 8 dpi for the historical variant, 10 dpi for the beta variant, at 12 dpi for variants delta and omicron (BA.1) and at 14 dpi for variant gamma (Fig 1 and S1 Table). In addition to viral RNA, infectious titers were measured for 41 animals. Infectious titers were measured by Tissue Culture Infectious Dose (TCID 50 ) from nasopharyngeal swab sampled at 2 time points per animal (day 2, 3 or 4 plus at day 5 or 7 post-infection). As we included several control animals from different studies, infected with either TCID 50 or Plaque Forming Units (PFU), all TCID 50 were converted to PFU assuming 1 PFU = 0.7 TCID 50 [START_REF] Davis | Nature of viruses[END_REF]. All infectious titers quickly dropped to undetectable levels for the historical variant at 5 dpi, where for the other variants the infectious titers remained consistent over the course of the infection (Fig 1). Fig 1 . Longitudinal measurements of genomic RNA, subgenomic RNA and infectious titers in 78 infected cynomolgus macaques. Both limit of quantification and detection are depicted as empty dots, the latter being lower. Upper limit of detection is depicted as filled squares.

Viral dynamic model

To account for the quick drop in infectious titers observed in the historical variant, (Fig 1 All models, except a model targeting the viral production parameter, provided an improvement of BIC compared to a target cell limited model (Table 1). We found that a model targeting the infectious ratio best described our data. In the following, we discuss the parameter values of the final constructed model accounting for both an effect of the immune effector and variant specific effect on the parameters (see below). For the historical variant, we estimated the infectivity rate parameter 𝛽 at 1.8610 -5 copies -1 .d -1 (95% confidence interval (CI) 110 -5 -3.3910 -5 ) and the loss rate of infected cells 𝛿 at 1.38 d -1 (95% CI 1.22 -1.55), corresponding to a half-life of 12 hours. We estimated the viral load production parameter 𝑝 at 9.4410 5 copies.cells -1 .day -1 (95% CI 2.110 5 -1.6810 6 ). This corresponds to a withinhost basic reproductive number 𝑅 0 (i.e., the number of newly infected cells by one infected cell at the beginning of the infection) of 3.1 (95% CI 2 -4.3) and a burst size (i.e the total number of infectious virus produced by one cell over its lifespan at the beginning of the infection) of 136 (95% CI 121 -153 ).

Model 4

Reduction of the viral production +9 Fig1. Schematic model of SARS-CoV-2 infection and action of the immune system. The basic model is a target cell limited model without any immune response. The parameters are : 𝜷 the infectivity rate, 𝒌 the transfer rate between non-productive and productive infected cells, 𝜹 the loss rate of productive infected cells, 𝒑 the viral production rate, 𝝁 the ratio of infectious virus, 𝒈 the transfer rate between the compartments of the immune response and 𝒄 the loss rate of both infectious and non-infectious virus

VoC specific effect on viral dynamic parameters

Once an effect of the immune response was selected, a covariate search algorithm was used to find the most likely VoC associated effects (see methods) and considered the historical variant as the reference. Several variant-specific covariates were found on viral kinetics parameters that we detail below (Fig 3 and S2 Table). First, beta variant was characterized with a reduced infected cells death rate (𝛿) by a factor of 0.7 (95% CI 0.6 -0.9) compared with the historical variant (p-value < 0.01). This led to an infected cell half-life of 17 hours and resulted in a longer period of viral load shedding as infected cells produced viruses for longer period of time. Gamma variant had an effect on the parameter 𝜃 (p-value < 0.001), the amount of immune effector 𝐹 20 required to reduce by half the infectious ratio, increasing it by a factor of 9508 (95% CI 387 -50 041) resulting in higher peak viral load and a longer duration of infectious virus shedding (Fig 4). Variant delta is characterized by an effect on both 𝜃 (p-value < 0.001) and the viral production parameter 𝑝 (p-value < 0.05), increasing those parameters by factors 336 (95% CI 49 -1191) and 1.78 (95% CI 1 -3) respectively. Finally, omicron variant (BA.1) affected the parameters of the immune system 𝜃 (p-value < 0.001), the viral production rate parameter 𝑝 (p-value < 0.001) and the infectious ratio 𝜇 (p-value < 0.001) modifying them by factors 229 (95% CI 27 -884), 0.07 (95% CI 0.02 -0.2) and 18 (95% CI 4 -51) respectively (Fig 4). The model well reproduced the viral load of all animals in the individuals fits (S2 Fig) . Additionally, we performed a sensitivity analysis on our best model (i.e. Model 1 including an effect on the infectious ratio μ). We tested several delays of the immune effector (from 1 to 6 days post infection) and several numbers of transfer compartments (from 5 to 30) and performed the covariate search on all models. We found that a delay of 3 days yielded the best results (S3 Table ) and very similar covariate were selected across all models (S3 and S4 Fig). We represent the mean value and 95% confidence interval of populations parameters for each variant. We represent only parameters having at least one variant-specific effect. Full table for population parameters is in S2 Using simulations, we sampled parameters considering both the uncertainty in the estimation and the inter-individual variability (see methods) .We represent the mean viral load of all variants and its 95% confidence interval. Dotted lines are the limits of detections

Predicted impact of variants in a natural infection setting

The main limitation of translating these results to humans is the fact that infection in animals is done with a large inoculum dose (10 5 -10 6 PFU), while human infections are presumably initiated with much lower virus dose [START_REF] Tang | Aerosol transmission of SARS-CoV-2? Evidence, prevention and control[END_REF]. Human experimental infections were performed with 10 TCID 50 [START_REF] Killingley | Safety, tolerability and viral kinetics during SARS-CoV-2 human challenge in young adults[END_REF] in the nose, i.e., 10,000-100,000 times less virus than in the animal model. Using simulations with lower inoculum, considering both uncertainty in the estimation and inter-individual variability (see methods), we are able to derive metrics of interest for each variant. The historical variant is characterized by a mean time to peak of 4.3 dpi (95% CI 3.7 -4.8) and of 3.5 dpi (95% CI 3 -3.9) for genomic RNA and infectious titers respectively. We found a mean peak viral load of 6.3 log 10 copies/mL (95% CI 5.5 -7) and of 2.1 PFU/mL (95% CI 1.2 -2.9) for genomic RNA and infectious titers, respectively. The reduced infected cell clearance rate of the beta variant resulted in a longer period of viral load shedding. The duration of the acute infection stage was consequently increased from 10.9 days (95% CI 9.5 -13.1) for the historical variant to 13.4 days (95% CI 11.1 -15.7) for the beta variant. All variants except beta have shown an effect on the antigen-mediated response, greatly reducing its impact on viral kinetics. As the effect of the antigen-mediated response was reduced, the infectious ratio was increased leading to more infectious particles produced over longer periods of time. This led to the increase of the infectious titers clearance stage duration from 1.5 days for the historical variant (95% CI 0.6 -1.9) to 6 days (95% CI 4.4 -7.5), 3.8 days (95% CI 3.1 -4.6) and 3.7 days (95% CI 2.8 -4.5) for the gamma, delta and omicron variants respectively (Fig 5). This is in line with numbers of studies showing the immune escape capabilities of those variants [START_REF] Chakraborty | A Detailed Overview of Immune Escape, Antibody Escape, Partial Vaccine Escape of SARS-CoV-2 and Their Emerging Variants With Escape Mutations[END_REF][START_REF] Mlcochova | SARS-CoV-2 B.1.617.2 Delta variant replication and immune evasion[END_REF][START_REF] Willett | SARS-CoV-2 Omicron is an immune escape variant with an altered cell entry pathway[END_REF]. An effect increasing the viral production parameter (𝑝), as observed for the delta variant, results in largely higher peak viral load of 7.6 log 10 copies/mL (95% CI 6.8 -8.2) and peak infectious titers of 3.7 PFU/mL (95% CI 2.4 -4.6). Conversely, an effect reducing the viral production parameter, as observed for the omicron variant, results in lower peak viral load compared to the historical variant of 5.6 log 10 copies/mL (4.8 -6.3) but very similar peak infectious titers at 2.8 PFU/mL (95% CI 1.9 -3.8). This is due to an effect of omicron on the infectious ratio, increasing the proportion of infectious virus produced. Discussion Here, we used mechanistic models to characterize in detail the viral dynamics of the main variants of concern in an experimental model of non-human primates. We evaluated the impact of an antigenmediated immune response on the viral dynamics and found that an effect reducing the infectious ratio best described our data. Some of the variants of concern, gamma, delta and omicron (BA.1) showed a strong ability to escape this response greatly increasing the number of infectious viruses produced over the course of the infection compared to the historical variant. Interestingly, the delta variant was associated with an increased viral production rate, whereas the omicron variant was associated with a lower viral production rate but a higher infectious ratio. Using simulations in a natural infection scenario, we found that omicron infections, relative to delta infections, are associated with lower peak viral RNA and reduced duration of viral RNA clearance while having similar peak infectious titers and duration of infectious titers clearance. These results suggest that omicron's infectiousness cannot be attributed to an increased viral RNA production but maybe due to an immune escape coupled with an increased infectious ratio, greatly increasing the number of infectious particles produced. Although many other factors are at play to explain the increased transmissibility of certain variants of concern, differences in viral dynamics can provides insights into the biology of those variants. As such, delta infections featuring increased peak viral load and infectious titers can increase the risk of "superspreading" events and infections outside of close-contact settings [START_REF] Kissler | Viral Dynamics of SARS-CoV-2 Variants in Vaccinated and Unvaccinated Persons[END_REF]. Omicron infections, on the other hand, featuring lower peak viral load concentration [START_REF] Sentis | SARS-CoV-2 Omicron Variant, Lineage BA.1, Is Associated with Lower Viral Load in Nasopharyngeal Samples Compared to Delta Variant[END_REF] but similar infectious titers respective to other variants, may result in transmission events that would not occur with other variants because insufficient infectious titers would be produced. These results are coherent with reports showing lower pathogenicity of omicron infection [START_REF]Severity of disease associated with Omicron variant as compared with Delta variant in hospitalized patients with suspected or confirmed SARS-CoV-2 infection[END_REF], as they are associated with lower viral burden. The combination of immune escape abilities, increased infectious ratio and longer duration of infectious virus shedding could be a possible mechanism to explain the enhanced transmissibility of omicron variant. As such, the quantification of infectious titers over time is crucial to inform further public health policies and adjust the isolation period accordingly. Our study has some important limitations. First, although we can characterise in detail the viral dynamics of SARS-CoV-2 in nonhuman primates in a controlled environment, the inoculated dose is extremely high (10 000 to 100 000 times higher ( 19)) compared with human infection. This leads to rapid saturation of target cells and makes it difficult to accurately estimate the early phase of infection. In the future, studies evaluating lower inoculum in NHP can greatly improve the precision in the estimation of the early phase of infection. Second, we developed an extension of the target cell limited model considering the effect of an antigen-mediated immune response decreasing the infectious ratio 𝜇. We here attribute this effect to the immune system but we have no information to which immune effectors (antibodies, cytokines, cytotoxic cells, natural killers, intracellular processes etc…) this could be linked if even attributable to one. This type of antigen-mediated response allows us to incorporate the effect of time on a parameter but the underlying biological mechanisms are unclear and may be due to inherent differences between variants not captured by any covariates. Third, we assumed a 3-days delay in the establishment of this antigen-mediated reduction of the infectious ratio and verified that it performed best in a sensitivity analysis (S4 Fig) . Although there is some variability, the covariates search is overall consistent. Fourth, the infectious titers are only a measure of in vitro infectivity, and to what extent they translate into infectiousness is unknown. In addition, both the upper and lower limit of quantification makes it difficult to precisely estimate the infectious ratio parameter 𝜇. Finally, in a context where more than half of the world population has received at least one dose of COVID-19 vaccine [START_REF] Ritchie | Coronavirus Pandemic (COVID-19)[END_REF], there is very little information on the natural infection with different variants. Additional data with vaccinated animals could help differentiate certain aspects of the abilities of the new variants to escape the immune system.

Materials and methods

Experimental procedure

Data comes from studies performed on cynomolgus macaques to evaluate the viral dynamics of SARS-CoV-2 variants. Our study includes 78 cynomolgus macaques (Macaca fascicularis) coming from control arms of several studies and have received no pharmacological interventions besides placebo. All animals were infected with doses ranging from 7×10 4 For each group both genomic RNA and subgenomic RNA swab samples were quantified using real time PCR in both the nasopharynx and in the trachea. For 41 animals (13 Historical, 3 Beta, 5 Gamma, 11 Delta and 7 Omicron (BA.1)) infectious titers were measured at 2 time points, early (2, 3 or 4 days post infection) and late (5 or 7 days post infection) using Tissue Culture Infectious Dose (TCID 50 ) from nasopharyngeal swab samples [START_REF] Antonio | SARS-CoV-2 viral dynamics in non-human primates[END_REF]. As we included animals from different studies that were inoculated with different methods (PFU or TCID 50 ), we normalized all measures of infectious titers by converting all TCID 50 measurements to Plaque Forming Units (PFU) using the formula 1 PFU = 0.7 TCID 50 [START_REF] Davis | Nature of viruses[END_REF]. As no infectious titers were measured in the trachea samples, we focused the main analysis on the nasopharyngeal compartment. The results mainly focus on the genomic viral load as the subgenomic is a directly proportional to the latter.

Basic viral dynamic model

We used a previously described model of SARS-COV-2 viral dynamics to reconstruct the nasopharyngeal viral load of infected animals. In this model, target cells (T) become infected cells (I 1 ) at a rate β. Infected cells transition into productive infected cells (I 2 ) at a rate k and produce infectious virus (V I ) at a rate 𝑝𝜇 and non-infectious virus (V NI ) at a rate 𝑝(1 -𝜇). Productive infected cells are cleared at a rate δ and both infectious and non-infectious virus are cleared at a rate c. The basic withinhost reproductive number, representing the number of newly infected cells by one infected cell, is 𝑅 0 = 𝛽𝑝𝑇 0 𝜇 𝑐𝛿 and the burst-size, representing the number of infectious virus produced by on infected cells over its lifespan, is 𝑁 = 𝑝𝜇 𝛿 . The model is described with the following set of ordinary differential equations:

𝑑𝑇 𝑑𝑡 = -𝛽𝑉 𝐼 𝑇 (1) 𝑑𝐼 1 𝑑𝑡 = 𝛽𝑉 𝐼 𝑇-𝑘𝐼 1 ( 2
)
𝑑𝐼 2 𝑑𝑡 = 𝑘𝐼 1 -𝛿𝐼 2 (3)
𝑑𝑉 𝐼 𝑑𝑡 = 𝑝𝜇𝐼 2 -𝑐𝑉 𝐼 (4)
𝑑𝑉 𝑁𝐼 𝑑𝑡 = 𝑝(1 -𝜇)𝐼 2 -𝑐𝑉 𝑁𝐼 (5)

Assumption on parameter values

Some parameters of the model were fixed to ensure identifiability. The transfer rate parameter between infected cells and productive infected cells was fixed to k = 4 day -1 (corresponding to a mean duration of the eclipse phase, i.e. the time for infected cells to start producing viruses, of 1 𝑘 = 6 hours) [START_REF] Ke | In vivo kinetics of SARS-CoV-2 infection and its relationship with a person's infectiousness[END_REF]. The viral clearance c was set to 10 day -1 based on previous work [START_REF] Maisonnasse | Hydroxychloroquine use against SARS-CoV-2 infection in non-human primates[END_REF][START_REF] Antonio | SARS-CoV-2 viral dynamics in non-human primates[END_REF][START_REF] Néant | Modeling SARS-CoV-2 viral kinetics and association with mortality in hospitalized patients from the French COVID cohort[END_REF]. As only the product 𝑝𝑇 0 is identifiable, we choose to fix the initial number of target cell to 𝑇 0 = 12 500 cells following the same assumptions as in [START_REF] Antonio | SARS-CoV-2 viral dynamics in non-human primates[END_REF]. As the nasal cavity of the animals is small, a substantial fraction of the inoculum does not penetrate the upper respiratory tract. To account for this, we introduced a parameter ℎ representing the proportion of the inoculum that arrive on the site of infection. We fixed this parameter at 20% with a standard deviation of 20% to allow for individual variability. As both the initial infectious inoculum and the number of RNA copies were known we used that information as our initial condition for the infectious virus and non-infectious virus compartment. Therefore, our initial conditions were set to :

𝑇 0 (𝑡 = 0) = 1.25 × 10 4 (6) 𝐼 1 (𝑡 = 0) = 0 (7) 𝐼 2 (𝑡 = 0) = 0 (8) 𝑉 𝐼 (𝑡 = 0) = 𝑉 𝐼 (0) 𝑖 × ℎ 𝑖 × 0.1 (9)
𝑉 𝑁𝐼 (𝑡 = 0) = (𝑉 𝑁𝐼 (0) 𝑖 -𝑉 𝐼 (0) 𝑖 ) × ℎ 𝑖 × 0.1 (10) Where 𝑉 𝐼 (0) 𝑖 is the administered dose in PFU of subject i, 𝑉 𝑁𝐼 (0) 𝑖 is the total number of RNA copies in the initial inoculum of subject i and h i is the proportion of the inoculum actively initiating the infection.

Models incorporating antigen-mediated immune response

To account for the quick drop in infectious titers observed for the historical variant (Fig 1 andS2 Fig), we tested several models incorporating an action of an antigen-mediated immune response. We assumed a delay of 3 days for the immune response to take place to account for the differentiation and proliferation of the immune response [START_REF] Janeway | Principles of innate and adaptive immunity[END_REF]. We modelled this delayed immune effector compartment using the Linear Chain Trick (LCT) assuming an Erlang distribution with 𝑗 = 20 transition compartment and a mean time spent in those compartment of 𝜏 = 3 d -1 [START_REF] Hurtado | Building mean field ODE models using the generalized linear chain trick & Markov chain theory[END_REF]. This number of compartments allowed us to shift the distribution of the time spent in the transition's states from an exponential to a normal distribution. The equations for the transfer compartments are written as follows: [START_REF] Hurtado | Building mean field ODE models using the generalized linear chain trick & Markov chain theory[END_REF] In the following only the compartment 𝐹 20 will serve as the effector for the action of the immune system. The transfer rate parameter 𝑔 is then written as 𝑗 𝜏 and fixed to 6.67 d -1 and the loss rate of the final effector 𝑑 𝐹 is fixed to 0.4 d -1 [START_REF] Néant | Modeling SARS-CoV-2 viral kinetics and association with mortality in hospitalized patients from the French COVID cohort[END_REF]. Several modes of action of the response system were tested: Model 1 : Immune effector decreases the infectious ratio 𝜇 In this model, the immune effector directly decreases the infectious ratio parameter 𝜇 using an Emax function type expression : (32) With 𝜃 being the amount of immune effector 𝐹 20 needed to reduce by half the infectious ratio. Model 2 : Immune effector increases infected productive cells death rate 𝛿 The death rate of infected cells is increased in proportion to the amount of immune effector 𝐹 20 .

𝑑𝐹 1 𝑑𝑡 = 𝐼 2 -𝑔𝐹 1 (11) 𝑑𝐹 2 𝑑𝑡 = 𝑔𝐹 1 -𝑔𝐹 2 (12) ⋮ 𝑑𝐹 20 𝑑𝑡 = 𝑔𝐹 19 -𝑑 𝐹 𝐹 20
𝑑𝑉 𝐼 𝑑𝑡 = 𝑝𝜇 ( 1 -𝐹 20 𝐹 20 + 𝜃 ) 𝐼 2 -𝑐𝑉 𝐼 (31) 𝑑𝑉 𝑁𝐼 𝑑𝑡 = 𝑝 ( 1 -𝜇 ( 1 -𝐹 20 𝐹 20 + 𝜃 ) ) 𝐼 2 -𝑐𝑉 𝑁𝐼

𝑑𝐼 2 𝑑𝑡

= 𝑘𝐼 1 -𝛿(1 + 𝜑𝐹 20 )𝐼 2 [START_REF] Comets | Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm[END_REF] Where 𝜑 is the strength of the immune system. Model 3 : Immune effector reduces the infectivity rate 𝛽 In this model, the immune effector blocks virus entry in the cells by reducing the infectivity parameter 𝛽. 𝑑𝑇 𝑑𝑡

= -𝛽(1 -𝜑𝐹 20 )𝑉 𝐼 𝑇 (34) 𝑑𝐼 1 𝑑𝑡 = 𝛽(1 -𝜑𝐹 20 )𝑉 𝐼 𝑇 -𝑘𝐼 1 (35)
Model 4 : Immune effector reduces the production rate 𝑝 In the same way as model 1, the viral load production parameter is reduced by the immune effector with an Emax type function:

𝑑𝑉 𝐼 𝑑𝑡 = 𝑝 ( 1 -𝐹 20 𝐹 20 + 𝜃 ) 𝜇𝐼 2 -𝑐𝑉 𝐼 (36) 𝑑𝑉 𝑁𝐼 𝑑𝑡 = 𝑝 ( 1 -𝐹 20 𝐹 20 + 𝜃 ) (1 -𝜇 )𝐼 2
-𝑐𝑉 𝑁𝐼 (37) All models were compared based on the Bayesian Information Criterion (BIC). We selected the model that yielded the lowest BIC and the best individual fits.

Statistical model

Parameter estimation was performed using non-linear mixed effect modelling. The statistical models describing the genomic RNA, subgenomic RNA and the infectious titers are:

𝑦 1 𝑖𝑗 = log 10 𝑉(𝑡 𝑖𝑗 ,Ψ 𝑖 ) + 𝑒 1 𝑖𝑗 (38) 𝑦 2 𝑖𝑗 = log 10 𝑓 × 𝐼 2 (𝑡 𝑖𝑗 , Ψ 𝑖 ) + 𝑒 2 𝑖𝑗 (39) 𝑦 3 𝑖𝑗 = log 10 𝑉 𝐼 (𝑡 𝑖𝑗 , Ψ 𝑖 ) + 𝑒 3 𝑖𝑗 (40)
Where the superscript 1, 2 and 3 refers to the genomic RNA, subgenomic RNA and infectious titers, respectively. We denote y ij is the j th observation of subject i at time t ij , with i ∈ 1, …, N and j ∈ 1, …, n i with N the number of subject and n i the number of observations for subject i. The function describing the total viral load kinetics 𝑉(𝑡 𝑖𝑗 , 𝛹 𝑖 ) predicted by the model at time t ij defined as: 𝑉 𝐼 (𝑡 𝑖𝑗 , 𝛹 𝑖 ) + 𝑉 𝑁𝐼 (𝑡 𝑖𝑗 , 𝛹 𝑖 ) predicted by the model at time t ij . The The vector of individual parameters of subject i is noted Ψ i and e ij is the additive residual Gaussian error of constant standard deviation σ. The vector of individual parameters depends on a fixed effects vector and on an individual random effects vector, which follows a normal centered distribution with diagonal variance-covariance matrix Ω. All parameters follow a lognormal distribution to ensure positivity except both parameters µ and h which follows logit-normal distribution and are bounded between 0 and 1. We assumed random effect on all parameters and removed them using backward procedure, if they were < 0.1 or their RSE > 50%. All biomarkers (i.e. genomic RNA, subgenomic RNA and infectious titers) were fitted simultaneously.

Selection of variant-specific effect on the viral dynamic parameters

Using the best model selected at the previous step, we sought to identify VoC-specific effect on the parameters of the model (𝛽, 𝛿, 𝑝,𝜇 and 𝜃). We first performed a backward selection of the random effects removing non-significant ones (i.e. relative standard error > 50%) if the BIC wasn't degraded by more than 2 points. We then used the Conditional Sampling use for Stepwise Approach on Correlation tests (COSSAC) to identify variant specific effect [START_REF] Ayral | A novel method based on unbiased correlations tests for covariate selection in nonlinear mixed effects models: The COSSAC approach[END_REF]. Then a backward procedure was used to remove any non-significant covariate effect with a Wald test (i.e. the covariate was removed if its coefficient effect relative standard error was > 50%). This procedure was repeated until all nonsignificant covariate effects had been eliminated. Additionally, we performed a sensitivity analysis on our best structural model. We tested for several delays in the establishment of the antigen-mediated effector (from 1 to 6 days) and on the number of transitions compartments (from 5 to 30) and then performed the covariate search on all model combinations .

Simulation of natural human infection

Finally, we used our final model to assess the impact of variants of concern on viral load and viral infectivity in a natural infection setting. We used a starting inoculum of 10 infectious virus, as described in an experimental challenge conducted in England [START_REF] Killingley | Safety, tolerability and viral kinetics during SARS-CoV-2 human challenge in young adults[END_REF] to simulate a human infection. The initial conditions are then written as: 𝑉 𝐼 (𝑡 = 0) = 10 𝑉 𝑁𝐼 (𝑡 = 0) = 0 We provided confidence interval on the mean predicted viral load, considering both the uncertainty in the estimation and the inter-individual variability. We first sampled M = 100 population parameters in their estimation distribution and then, for each variant, sampled N = 30 individual parameters from each sets of population parameters (leading to 3000 individual parameters per variant). We calculated the predicted viral load of all individuals and derived the mean viral load over the simulated individuals at all times with its 95% inter quantile range. Additionally, we provided the distribution of several viral dynamic metrics, namely:

-the area under viral load curve, -the peak and time to peak viral load -the duration of the clearance stage, calculated as the time interval between the peak viral load and the time to undetectable viral load -the duration of the acute phase, calculated as the time between the first and the last detectable viral load [START_REF] Kissler | Viral dynamics of acute SARS-CoV-2 infection and applications to diagnostic and public health strategies[END_REF].

Parameter estimation

All parameters were estimated by computing the maximum-likelihood estimator using the stochastic approximation expectation-maximization (SAEM) algorithm implemented in Monolix Software 2020R1 [START_REF] Comets | Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm[END_REF][START_REF]Lixoft[END_REF]. Standard errors and the likelihood were computed by importance sampling.

Acknowledgement

We would like to thank everyone in the CEA and at Pasteur Institute that have helped for data collection. We thank Alan Perelson for helpful discussions.

Figures Fig 1.

Longitudinal measurements of genomic RNA, subgenomic RNA and infectious titers in 78 infected cynomolgus macaques. Both limit of quantification and detection are depicted as empty dots, the latter being lower. Upper limit of detection is depicted as filled squares (present only un infectious titers). We represent the mean value and 95% confidence interval of populations parameters for each variant. We represent only parameters having at least one variant-specific effect. Full table for population parameters is in S2

  and S1 Fig) several models incorporating an action of an antigen-mediated immune response were tested (Fig 2).

Fig 3 .

 3 Fig 3. Estimated population parameters for each variant.We represent the mean value and 95% confidence interval of populations parameters for each variant. We represent only parameters having at least one variant-specific effect. Full table for population parameters is in S2 Table. The dashed black line represents the historical value.

Fig 4 .

 4 Fig 4. Simulation of variant of concern impact on viral load.Using simulations, we sampled parameters considering both the uncertainty in the estimation and the inter-individual variability (see methods) .We represent the mean viral load of all variants and its 95% confidence interval. Dotted lines are the limits of detections

Fig 5 .

 5 Fig 5. Impact of VoC on viral load metrics in the context of an infection with a low inoculum.We represent the mean and 95% confidence interval for each variant. The dashed black line represents the historical mean value.

  to 10 6 PFU of different SARS-CoV-2 strains. Animals are infected via both nasopharyngeal and intratracheal route with 10% of the initial volume administered in the nose and 90% in the trachea. The study is composed of 5 groups, each infected with a different SARS-CoV-2 strains: 44 Historical (hCoV-19/France/lDF0372/2020 strain; GISAID EpiCoV platform under accession number EPI_ISL_406596), 9 Bêta (B.1.351 -hCoV-19/USA/MD-HP01542/2021, BEI NR-55283), 5 Gamma (P.1 -hCoV-19/Japan/TY7-503/2021, BEI NR-54984), 11 Delta (B.1.617.2 -hCoV-19/USA/MD-HP05647/2021, BEI NR-55674) and 9 Omicron (B.1.1.529 -hCoV-19/USA/MD-HP20874/2021, BEI NR-56462).

Fig2.Fig 3 .

 3 Fig2. Schematic model of SARS-CoV-2 infection and action of the immune system. The basic model is a target cell limited model without any immune response. The parameters are : 𝜷 the infectivity rate, 𝒌 the transfer rate between non-productive and productive infected cells, 𝜹 the loss rate of productive infected cells, 𝒑 the viral production rate, 𝝁 the ratio of infectious virus, 𝒈 the transfer rate between the compartments of the immune response and 𝒄 the loss rate of both infectious and non-infectious virus.

Fig 4 .θ 5 .

 45 Fig 4. Simulation of variant of concern impact on viral load.Using simulations, we sampled parameters considering both the uncertainty in the estimation and the inter-individual variability (see methods) .We represent the mean viral load of all variants and its 95% confidence interval. Dotted lines are the limits of detections.

Table 1 : Alternative immune response models.

 1

	Models	Description	ΔBIC
	Reference model	Absence of immune response	-
	Model 1	Reduction of the infectious ratio	-𝟒𝟐.𝟖
	Model 2	Increase in infected cell clearance	-14.8
	Model 3	Reduction of viral infection rate	-36.1

  Table. The dashed black line represents the historical value.

  Table. The dashed black line represents the historical value.

Table. Characteristics of the 78 animals analysed

  . Descriptive statistics of the animals calculated on the raw data. S2

	S1
										Genomic RNA
	Area under viral load curve		Time to peak		Peak viral load		Clearance stage		Acute stage
		(log10 copies.days/mL)	viral load (days)		(log10 copies/mL)		duration (days)		duration (days)
	Omicron
	Delta
	Gamma
	Beta
	Historical
	70	80	90 100 110	6	9	12	15	5	6	7	8	6	8	10	10	12	14	16	18
										Infectious titers
	Area under viral load curve		Time to peak		Peak viral load		Clearance stage		Acute stage
		(log10 copies.days/mL)	viral load (days)		(log10 copies/mL)		duration (days)		duration (days)
	Omicron
	Delta
	Gamma
	Beta
	Historical
		25	30	35	40	5	10		15	2	3	4		2	4	6	2.5 5.0 7.5 10.0 12.5

Table. Estimates of the population parameter and covariate effects for the best model

  . **The standard error for the 𝑅 0 parameters were calculated using the delta method.

	461
	462	S3 Table.
		Strains	Number of animals	Mean weight (kg)	Mean peak viral load (log 10 copies/mL)	Mean peak PFU (log 10 PFU/mL)	Mean time to first undetectable viral load	Mean time to first undetectable PFU
		Historical	44	3.7		7.6	2.3	8	4
		Beta	9	4.9		7.1	3.2	10	6
		Gamma	5	4.2		7.8	3	14	5
		Delta	11	3.6		8.1	2.9	12	5
		Omicron	9	4.6		6.4	2.4	12	7
		Population parameters (unit)		Fixed effect (RSE%)	SD of random effect (RSE%)
		𝛽 (copie.d -1 )			1.85 × 10 -5 (33)	-
		𝑝 (copies.cell -1 .d -1 )			9.44 × 10 5 (40)	0.61 (17)
		𝛿 (d -1 )				1.38 (6)		0.2 (20)
		𝑓 (unitless)			1.36 × 10 -3 (19)	-
		𝜇 (unitless)			1.98 × 10 -4 (47)	-
		𝜃 (unitless)			0.19 (45)		0.32 (103)
		Covariate model			Covariate effect (RSE%)	p-value
		Beta on 𝛿			-0.357 (32)		0.00201
		Gamma on 𝜃			8.39 (15)		< 10 -6
		Delta on 𝑝			0.554 (50)		0.047
		Delta on 𝜃			5.49 (15)		< 10 -6
		Omicron on 𝑝			-2.82 (19)		< 10 -6
		Omicron on 𝜇			2.7 (23)		1.98 × 10 -5
		Omicron on 𝜃			5.04 (18)		< 10 -6
		Basic reproductive number			Value (RSE%)**
		𝑅 0				3.1 (19)
		𝑅 0𝑏𝑒𝑡𝑎				4.5 (20)
		𝑅 0𝑔𝑎𝑚𝑚𝑎			3.1 (19)
		𝑅 0𝑑𝑒𝑙𝑡𝑎			5.4 (34)
		𝑅 0𝑜𝑚𝑖𝑐𝑟𝑜𝑛			2.8 (28)
		Residual errors			Value (RSE%)
		𝜎 𝐺𝑒𝑛𝑜𝑚𝑖𝑐 𝑅𝑁𝐴 (log 10 copies/mL)		0.98 (4)
		𝜎 𝑆𝑢𝑏𝑔𝑒𝑛𝑜𝑚𝑖𝑐 𝑅𝑁𝐴 (log 10 copies/mL)	0.89 (6)
		𝜎 𝐼𝑛𝑓𝑒𝑐𝑡𝑖𝑜𝑢𝑠 𝑡𝑖𝑡𝑒𝑟𝑠 (log 10 PFU/mL)		1.79 (14)

Sensitivity analysis on the delayed immune response

  . Using the best structural model (i.e. 463Model 1 including an effect on the infectious ratio) we tested several delays for the immune response to 464 take place and performed the covariate search algorithm on all models.

	Model ID	Number of transfer compartments 𝐹	Delay (days)	Transfer rate parameter 𝑔 (d -1 )	BIC before COSSAC	BIC after COSSAC (ΔBIC )
	1	5	1	5	2451	(-22)
	2	5	2	2.5	2409	(-25)
	3	5	3	1.666666667	2405	(-31)
	4	5	4	1.25	2408	(-35)
	5	5	5	1	2408	(-33)
	6	5	6	0.833333333	2409	(-30)
	7	10	1	10	2432	(-23)
	8	10	2	5	2409	(-36)
	9	10	3	3.333333333	2410	(-49)
	10	10	4	2.5	2411	(-30)
	11	10	5	2	2413	(-47)
	12	10	6	1.666666667	2414	(-47)
	13	20	1	20	2426	(-24)
	14	20	2	10	2409	(-46)
	15	20	3	6.666666667	2411	(-51)
	16	20	4	5	2414	(-37)
	17	20	5	4	2416	(-37)
	18	20	6	3.333333333	2417	(-36)
	19	30	1	30	2424	(-27)
	20	30	2	15	2408	(-31)
	21	30	3	10	2413	(-28)
	22	30	4	7.5	2417	(-54)
	23	30	5	6	2419	(-39)
	24	30	6	5	2420	(-27)
	465
	466
	467
	468
	469
	470
	471
	472

Supplementary material

We performed a sensitivity analysis on our best model. The model IDs are represented on top, as described in S3 Table . The scale represents the magnitude of the covariate effect rescaled for each row with 0 being the minimum value and 1 the maximum. Empty tiles indicate that no covariates were selected for this variant-parameter relationship. S4 Fig. Consistency of the covariate selection algorithm. We represent the number of times a covariate was found on a variant-parameter relationship across all 24 models. Empty tiles indicate that no covariates were found for this variant-parameter relationship.


Marjan Boerma

; Omid Azimzadeh

Nadia Pasinetti

Virginie Monceau

Editorial : "Cardiotoxicity Induced by Radiotherapy and/or Chemotherapy After Cancer Treatment "

Keywords: advanced, targeted therapy, cardiotoxicity, radiotherapy, cancer treatment

It is our privilege to present 11 articles in this Frontiers Oncology Research Topic on "cardiotoxicity induced by radiotherapy and/or chemotherapy after cancer treatment". The therapeutic management of cancer has progressed over the last 10 years and has led to a significant increase in patient survival rates. Radiotherapy (RT) is classically applied in multiple fractions administered over several weeks to kill cancer while sparing normal tissue as much as possible. However, the corollary is the development of delayed toxicities that alter the quality of life of patients. Chemo-and/or RT-induced cardiovascular disease (CVD) is recognized as a worrisome side effect in patients with thoracic cancers. This toxicity appears to be amplified if there are comorbidities as presented in the article by Kathryn Banfill et al. [START_REF] Banfill | Heart dose and cardiac comorbidities influence death with a cardiac cause following hypofractionated radiotherapy for lung cancer[END_REF]. that shows the effect of radiation dose on cardiac death after thoracic RT is different in patients with and without cardiac comorbidities. Therefore, CVD risk factors should be identified and managed in conjunction with RT for lung cancer. In the same context, Ciprian Camil Mirestean et al. describe the concerns that should be considered in breast cancer (BC) therapy. BC is the most common cancer in women worldwide, often treated with RT including whole breast irradiation (WBI). Recent data recommend the use of hypofractionation (HF)-WBI regardless of age or stage of disease. However, some studies reported an increased incidence of cardiac death with HF-WBI, particularly in patients with pre-existing cardiac risk factors at the time of treatment. There is also a need to develop multivariable radiobiologic models including histologic, molecular, and clinical parameters to identify risk groups and dosimetric tolerance to limit the incidence of late cardiac events [START_REF] Mirestean | Hypofractionated Whole-Breast Irradiation Focus on Coronary Arteries and Cardiac Toxicity-A Narrative Review[END_REF]. Azadeh Abravan et al. describe the use of real-world data (RWD) to study radiation-induced heart disease (RIHD) in lung cancer patients, summarizing how heart dosimetric factors associate with outcome, strength, and limitations of the RWD studies, and how RWD can be used to assess a change to cardiac dose constraints [START_REF] Abravan | Role of Real-World Data in Assessing Cardiac Toxicity After Lung Cancer Radiotherapy[END_REF]. Since RT of BC can result in an increased risk of long-term major CVD events, it is critical to detect subclinical left ventricular (LV) dysfunction early in RT treated BC patients to determine the dose-response relationships between cardiac doses and these events. The results of Médéa Locquet et al. highlight that all cardiac doses were strongly associated with the development of subclinical LV dysfunction 6 months after RT. It remains to be determined whether global longitudinal strain measurements at baseline and 6 months after RT can predict subclinical events occurring 24 months after RT. In addition, few studies suggested that RT of BC can induce cardiac arrhythmias and conduction disorders [START_REF] Locquet | Subclinical Left Ventricular Dysfunction Detected by Speckle-Tracking Echocardiography in Breast Cancer Patients Treated With Radiation Therapy: A Six-Month Follow-Up Analysis (MEDIRAD EARLY-HEART study)[END_REF]. However, the association between mean heart dose and doses to cardiac substructures is less studied. Mohamed Yassir Errahmani et al. performed an exploratory investigation on BC patients treated with RT, which is the first study suggesting that irradiation of the right atria may require special attention regarding the risk of cardiac arrhythmia and conduction disorders [START_REF] Errahmani | Association Between Cardiac Radiation Exposure and the Risk of Arrhythmia in Breast Cancer Patients Treated With Radiotherapy: A Case-Control Study[END_REF]. Regardless of the technique used, during RT, healthy tissues located in the irradiation fields are exposed. This is the main reason for the development of RT techniques in the treatment of left-sided BC, where conventional RT leads to cardiac radiation exposure, such as the deep inspiration breathhold (DIBH), as described by Yongkai Lu et al. Whether the DIBH regimen is an optimal solution for left-sided BC remains unclear. This meta-analysis aimed to elucidate the differences between DIBH and free-breathing for patients receiving RT for left-sided BC and provide a practical reference for clinical practice. This study suggests more widely use of DIBH in clinical practice because of its excellent dosimetric performance [START_REF] Lu | Comparison of Deep Inspiration Breath Hold Versus Free Breathing in Radiotherapy for Left Sided Breast Cancer[END_REF]. An article by Eva Kimpe et al. deals with the social impact of long-term effects of radiation-induced cardiotoxicity as an important concern during the treatment of BC . This study applied health economic modelling techniques to estimate attributed CVD-related costs and disutility. Their analyses suggest that the cost of past investments to achieve the mean heart dose (MHD) in current practice seems justified when considering the gains from cost resulting from radiation-induced CVD events [START_REF] Kimpe | Quantifying Societal Burden of Radiation-Induced Cardiovascular Events in Breast Cancer Survivors[END_REF]. The question remains to be answered is whether costs associated with further investments in technological advancements offset the expected benefit from reducing the MHD. Altogether, epidemiological and clinical data underline the importance of cardiac side effects after RT, but the pathophysiological, cellular and molecular bases of such side effects remain poorly understood. This Research Topic highlights research that creates new knowledge of cellular and molecular signaling pathways using preclinical animal models. To address the impact of genetic bases and disparate responses to RT, Saryleine Ortiz De Choudens et al. used Dahl salt-sensitive rats as a surrogate hypertension model to analyse the injury to the heart by irradiation [START_REF] Ortiz De Choudens | Lisinopril Mitigates Radiation-Induced Mitochondrial Defects in Rat Heart and Blood Cells[END_REF]. Similarly using a rat model of lung irradiation, Julia Wiedemann et al. demonstrated a series of pathologies following remodeling of the pulmonary vasculature. The similarities between the mechanisms of vascular remodeling in pulmonary hypertension and those after irradiation could be translated into interventions that benefit patients treated for thoracic tumors, where radiation to lung tissue often cannot be avoided [START_REF] Wiedemann | Radiation-induced cardiac side-effects: The lung as target for interacting damage and intervention[END_REF]. This knowledge greatly facilitates the discovery of biomarkers for cardiovascular toxicity, the identification of new cardioprotective therapeutic targets, and the optimization of prevention and intervention strategies in chemo-and RT. Cancer immune checkpoint inhibitors have led to recent advances in the field of cancer immunotherapy improving overall survival in multiple cancers with historically abysmal prognoses. Cardiac-specific immune-related adverse events are potentially fatal. However, the understanding of autoimmune cardiotoxicity remains limited due to its rareness. Omoruyi Credit Irabor et al. provide a literature review on the pathologic mechanisms, diagnosis, and management of autoimmune cardiotoxicity resulting from ICIs and offer a perspective on potential strategies and research to prevent and mitigate their occurrence [START_REF] Irabor | Overcoming the cardiac toxicities of cancer therapy immune checkpoint inhibitors[END_REF]. Thoracic RT has been associated with increased cardiac morbidity and mortality in numerous studies, including the landmark Darby´s study (2013), demonstrating a linear increase in cardiac mortality with increasing mean cardiac radiation dose. However, the extent to which cardiotoxicity has been incorporated as an endpoint in RT studies has not been standardized and is sometimes unknown. To better characterize cardiac toxicities, future prospective studies involving thoracic RT should include cardiotoxicity endpoints as recommended in the study by Rahul Neal Prasad et al. [START_REF] Prasad | Lack of Cardiotoxicity Endpoints in Prospective Trials Involving Chest Radiation Therapy: A Review of Registered, Latter-Phase Studies[END_REF]. We hope that this Research Topic will arouse the interest of radiation researchers, epidemiologists and clinicians to continue to pursue research that increases our knowledge on CVD following chemo-and RT, eventually implementing practices that will improve the safety of cancer therapy.

Authors Contribution: All authors listed have contributed equally. VM drafted the manuscript and incorporated the ideas of all authors. MB, OA and NP provided comments and approved of the final version.

Acknowledgments:

The authors would like to thank all scientists who contributed to this Research Topic as well as various reviewers/editors of the respective manuscripts, for their efforts, timely responses and enthusiasm. We also thank the Frontiers Editorial Office for their assistance and support.


Example of neighborhood : bit strings

Sample the neighborhood to have information on local features of the search space

From local information : deduce some global features like general shape of search space, "difficulty", etc.

S. Verel

Fitness landscapes and graphs

Goal of the fitness landscapes study

Study of the geometry of the landscape allows to study the difficulty, and design a good optimisation algorithm Fitness landscape is a graph (S, N , f ) where the nodes have a value (fitness) : can be "pictured" as a "real" landscape

Multimodal Fitness landscapes

Local optima s * :

no neighbor solution with higher fitness value

∀s ∈ N (s * ), f (s) < f (s * )
Search space

Fitness

S. Verel Fitness landscapes and graphs

Multimodal Fitness landscapes

Adaptive walk : (s 0 , s 1 , . . .) where s i+1 ∈ N (s i ) and f

(s i ) < f (s i+1 ) Hill-Climbing (HC) algorithm Choose initial solution s ∈ S repeat choose s ∈ N (s) such that f (s ) = max x∈N (s) f (x) if f (s) < f (s ) then s ← s end if until s is a Local optimum Basin of attraction of s * : {s ∈ S | HillClimbing (s) = s * }.

S. Verel

Fitness landscapes and graphs

Multimodal Fitness landscapes

Search space

Fitness

Optimisation difficulty : number and size of attractive basins (Garnier et al [10])

The idea :

if the size of attractive basin of global optima is relatively "small" the problem is difficult to optimize

The measure :

Length of adaptive walks (distribution, avg, etc.)

ρ(n) = E [(f (s i ) -f )(f (s i+n ) -f )] var (f (s i )) autocorrelation length τ = 1 ρ(1)

Neutrality and difficulty

We know for certain that :

No information is better than Bad information : Hard trap functions are more difficult than needle-in-a-haystack functions

Good information is better than No information

When there is No information : you should have a good method to find it ! S. Verel Fitness landscapes and graphs

In the following Neutral networks (Schuster 1994 [25])

genotypes space

Fitness

Fitness

Neutral Network Doors S. Verel Fitness landscapes and graphs

Definitions

Test of neutrality

isNeutral : S × S → {true, false} For example, isNeutral(s 1 , s 2 ) is true if : f (s 1 ) = f (s 2 ). |f (s 1 ) -f (s 2 )| ≤ 1/M with M is the search population size. |f (s 1 ) -f (s 2 )| is under the evaluation error.
Neutral neighborhood of s is the set of neighbors which have the same fitness f (s)

N neut (s) = {s ∈ N (s) | isNeutral(s, s )}

Neutral degree of s

Number of neutral neighbors : nDeg (s) = (N neut (s) -{s}).

S. Verel Fitness landscapes and graphs

Definitions

Neutral walk

W neut = (s 0 , s 1 , . . . , s m ) for all i ∈ [0, m -1], s i+1 ∈ N (s i ) for all (i, j) ∈ [0, m] 2 , isNeutral(s i , s j ) is true. Neutral Network graph G = (N, E )
N ⊂ S : for all s and s from V , there is a neutral walk belonging to V from s to s ,

(s 1 , s 2 ) ∈ E if they are neutral neighbors : s 2 ∈ N neut (s 1 )
A fitness landscape is neutral if there are many solutions with high neutral degree.

Autocorrelation of evolvability :

Evolvability evol = avg fitness in the neighborhood Autocorrelation of (evol(s 0 ), evol(s 1 ), . . .).

Informations :

if

Basins of attraction in combinatorial optimisation

Example of small NK landscape with N = 6 and K = 2

Color represent fitness value

Basins of attraction in combinatorial optimisation

Example of small NK landscape with N = 6 and K = 2

Color represent fitness value

Basins of attraction in combinatorial optimisation

Example of small NK landscape with N = 6 and K = 2

Each color corresponds to one basin of attraction

Basins of attraction are interlinked and overlapped

Basins have no "interior" S. Verel Fitness landscapes and graphs

Basins of attraction in combinatorial optimisation

Example of small NK landscape with N = 6 and K = 2

Basins of attraction are interlinked and overlapped ! Most neighbours of a given solution are outside its basin

Basin of attraction

Hill-Climbing (HC) algorithm

Choose initial solution s ∈ S repeat choose s ∈ N (s) such that f (s ) = max x∈N (s) f (x) if f (s) < f (s ) then s ← s end if until s is a Local optimum
Basin of attraction of s * : {s ∈ S | HillClimbing (s) = s * }.

S. Verel

Fitness landscapes and graphs

Methods

Extracted and analysed networks

(i) = 1 s i (k i -1) j,h w ij + w ih 2 a ij a jh a hi
where s i = j =i w ij , a nm = 1 if w nm > 0, a nm = 0 if w nm = 0 and k i = j =i a ij .

Disparity dishomogeneity of nodes with a given degree

Y 2 (i) = j =i w ij s i 2 S.

  Search space : S = {0, 1} N Algorithm : simple GA, hill-climbing, or simulated annealing, etc. {y ∈ S | d Hamming (y , x) ≤ 1} S. Verel Fitness landscapes and graphs Example of neighborhood : permutations Traveling Salesman Problem : find the shortest tour which cross one time every town Search space : S = {σ | σ permutations } Algorithm : simple EA operator : 2-opt N (x) = {y ∈ S | IP(y = op 2opt (x)) > 0} S. Verel Fitness landscapes and graphs Example of neighborhood Traveling Salesman Problem : find the shortest tour which cross one time every town Search space : S = {σ | σ permutations } Algorithm : simple EA operators : 2-opt and 3-opt N (x) = {y ∈ S | IP(y = op 2opt (x)) > 0 or IP(y = op 3opt (x)) > 0} S. Verel Fitness landscapes and graphs Example of neighborhood : memetic algorithms Algorithm : memetic algorithm, EA + operator hill-climbing N (x) = {y ∈ S | y = op HC (x)} Algorithm : memetic algorithm, EA + operator hill-climbing and bit-flip mutation 2 possibilities : Study 2 landscapes : one for HC operator, one for bit-flip mutation Study 1 landscape : N (x) = {y ∈ S | y = op HC (x) or IP(y = op bit-flip (x)) > } It depends on what you want to know S. Verel Fitness landscapes and graphs Point of view : Before putting a particular heuristic FL = (Sol., Neighbors, Fitness)

  F, IP) : probability space op : S → S stochastic operator of the local search X (s) = f (s) Y (s) = f (op(s)) Fitness Cloud of op Conditional probability density function of Y given X S. Verel Fitness landscapes and graphs Fitness cloud : Measure of evolvability Fitness f(s) Fitness f(op(s)) Neutral sets (done) : set of solutions with the same fitness ⇒ No structure Fitness cloud (done) : Bivariate density (f (s), f (op(s))) ⇒ Neighborhood relation between neutral sets Neutral networks (to be done) : ⇒ Neighborhood structure into the neutral sets : Graph S. Verel Fitness landscapes and graphs

Neighborhood

  nodes are connected if the energy barrier separating them is sufficiently low (transition state) S. Verel Fitness landscapes and graphs Basins of attraction in combinatorial optimisation Example of small NK landscape with N = 6 and K = 2 Bit strings of length N = 6 2 6 = 64 solutions one point = one solution S. Verel Fitness landscapes and graphsBasins of attraction in combinatorial optimisationExample of small NK landscape with N = 6 and K = 2Bit strings of length N = 6

  , 470 34 162139 4592 0.53 0.0070 0.050 0.0006 12 2, 254 32 227912 2670 0.44 0.0031 0.043 0.0003 14 3, 264 29 290732 2056 0.38 0.0022 0.040 0.0003 15 3, 868 33 321203 2061 0.35 0.0022 0.039 0.0004 Clustering Coefficient : For high K, transition between a given pair of neighboring basins is less likely to occur Disparity : For high K the transitions to other basins tend to become equally likely, an indication of the randomness of the landscape Adami, R. K. Belew, H. Kitano, and C. Taylor, editors, ALIFE VI, Proceedings of the Sixth International Conference on Artificial Life, pages 18-27. ALIFE, The MIT Press, 1998. Proceedings of the 2001 Congress on Evolutionary Computation CEC2001, pages 30-37, COEX, World Trade Center, 159 Samseong-dong, Gangnam-gu, Seoul, Korea, 27-30 2001. IEEE Press. U. Bastolla, M. Porto, H. E. Roman, and M. Vendruscolo. Statiscal properties of neutral evolution. Journal Molecular Evolution, 57(S) :103-119, August 2003. The roles of mutation, inbreeding, crossbreeding, and selection in evolution. In Proceedings of the Sixth International Congress of Genetics 1, pages 356-366, 1932.

	General network statistics N = 16 Weights distribution : transition probability between basins Interior and border size Shortest path length between local optima Shortest path length to global optima L. Barnett.
	K Ruggedness and neutrality -the NKp family of fitness nodes edges C w Ȳ landscapes.	d
	Verel 9129 2930 516 358 41791 4690 0.79 0.0154 0.084 0.0028 Fitness landscapes and graphs 0.96 0.0245 0.326 0.0579 56 14 0.92 0.0171 0.137 0.0111 93384 4394 0.65 0.0102 0.062 0.0011 Fitness landscapes and graphs S. Verel Fitness landscapes and graphs S. Verel Fitness landscapes and graphs S. Verel Fitness landscapes and graphs S. Verel Fitness landscapes and graphs 10 1S. Verel 2 33 15 4 178 33 6 460 29 8 890 33 0.0001 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.001 0.01 0.1 1 W K=2 K=4 K=10 K=14 K=17 distribution of the network P(wij=W) weights w ij for outgoing edges border with j = i in log-x scale, N = 18 Weights are small For high K the decay is faster Low K has longer tails On average, the transition probabilities are higher for low K (less local optima) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 2 4 6 8 10 12 14 16 18 average of the mean size K N=14 N=16 N=18 Average of the mean size of basins interiors Question : Do basins look like a "montain" with interior and border ? solution is in the interior if all neighbors are in the same basin Interior is very small Nearly all solution are in the K Average distance (shortest path) between nodes K up to K = 10, then stagnates optimum from all the other basins increasing K Average path length to the global Increase steadily with For a given N, increase with K Answer 0 50 100 150 200 250 300 2 4 6 8 10 12 14 16 18 average path length N=14 N=16 N=18 Question : Are the basins "far" from each other ? Increase with N ( of nodes 50 100 150 200 250 Question : Is the global optimum basin is far ? In C. Lionel Barnett. Netcrawling -optimal evolutionary search with neutral average path length to the optimum N=14 N=16 networks. More relevant for increases exponentially) 0 2 4 6 8 10 12 14 16 18 N=18 optimisation In

S. Verel

Fitness landscapes and graphs

T. Jones. Evolutionary Algorithms, Fitness Landscapes and Search. PhD thesis, University of New Mexico, Albuquerque, 1995. S. A. Kauffman. The Origins of Order. Oxford University Press, New York, 1993. M. Kimura. The Neutral Theory of Molecular Evolution. Cambridge University Press, Cambridge, UK, 1983. J. Lobo, J. H. Miller, and W. Fontana. Neutrality in technology landscape, 2004. M. Newman and R. Engelhardt. Effect of neutral selection on the evolution of molecular species. S. Verel Fitness landscapes and graphs

Overview and Motivation

Bring the tools of complex networks analysis to the study the structure of combinatorial fitness landscapes Goals : Understand problem difficulty, design effective heuristic search algorithms Methodology : Extract a network that represents the landscape (Inspiration from energy landscapes [START_REF] Verel | The network topology of a potential energy landscape : a static scale-free network[END_REF]

Weights of edges

From each s and s , p(s → s ) = IP(s = op(s))

For example, S = {0, 1} N and bit-flip operator

w ij : Total probability of going from basin b i to basin b j is the average over all s ∈ b i of the transition prob. to s ∈ b j :