R Gagnieu

Rodolphe Garraffo

C Guilhaumou

Solas

Nicolas Venisse
email: nicolas.venisse@chu-poitiers.fr

Gilles Peytavin

Stephane Bouchet

Marie-Claude Gagnieu

Romain Guilhaumou

Caroline Solas
email: caroline.solas@ap-hm.fr

Concerns about pharmacokinetic (PK) and pharmacokinetic-pharmacodynamic (PK-PD) studies

Keywords: COVID-19, Pharmacokinetics, Pharmacodynamics, PK-PD 3

Introduction

In the context of the COVID-19 pandemic, several existing approved drugs and experimental antiviral agents have been repurposed as potential antiviral candidates for the treatment of COVID-19 infection. While preliminary choices were essentially based on in vitro potency, clinical translation into effective therapies may be challenging due to unfavorable in vivo pharmacokinetic (PK) properties (i.e. plasma protein binding, tissue distribution, drug interactions) at the doses chosen for this new indication of COVID-19 infection. The particular conditions of COVID-19 infection (cytokine storm, multi-visceral failure and lifethreatening prognosis), patient co-morbidities (i.e. obesity, diabetes, cardiovascular complications) and the requirement for a short and rapidly effective treatment further complicate the choice of the ideal candidate.

Remdesivir, chloroquine derivatives (essentially hydroxychloroquine (HCQ) due to a better safety profile than chloroquine) and the anti-HIV agent lopinavir (LPV) were among the first to be tested due to an in vitro antiviral activity demonstrated against SARS-CoV-2 or other similar respiratory viruses (i.e. SARS-CoV, MERS-CoV). Although remdesivir is not yet commercially available, other agents are already easily accessible to clinical investigators as they are part of the therapeutic armamentarium of other diseases (i.e. systemic lupus erythematosus (SLE) for HCQ and human immunodeficiency virus (HIV) treatment for LPV) probably explaining the large number of ongoing clinical trials worldwide. Favipiravir, ribavirin, tocilizumab, ivermectin, nafamostat and other agents have also been proposed for treatment of COVID-19 infection either as antivirals or immunomodulatory agents [START_REF] Sanders | Pharmacologic Treatments for Coronavirus Disease 2019 (COVID-19): A Review[END_REF].

In the area of infectious diseases and antiviral drugs, pharmacological properties are of particular importance for treatment choices, evaluation and optimization. Indeed, suboptimal antiviral response may be a consequence of inadequate exposure and/or poor PK-PD properties of the studied drug. From the HIV pandemic, we have learned that maintaining sufficient plasma drug exposure is critical to stop virus replication and avoid emergence of resistances [START_REF] González De Requena | Nevirapine plasma exposure affects both durability of viral suppression and selection of nevirapine primary resistance mutations in a clinical setting[END_REF]. This has led to the implementation of strategies to optimize dosing regimen such as the PK "boosting", used in the lopinavir/ritonavir (LPV/r) association. In this area, increasing knowledge on the PK and pharmacokinetic-pharmacodynamic (PK-PD) relationships of antiretrovirals has also demonstrated its usefulness in treatment optimization through the use of therapeutic drug monitoring [START_REF] Boffito | Current status and future prospects of therapeutic drug monitoring and applied clinical pharmacology in antiretroviral therapy[END_REF]. In the context of COVID-19 infection, optimizing drug exposure at the site of infection, i.e. in the respiratory tract, is probably the key to successful treatment.

Accurate collection of PK and PD data is therefore of primary importance, especially for these repurposed drugs. We believe that extrapolation of PK data from other clinical situations may require specific caution due to different physiopathological conditions. In the context of the current global emergency, the number of clinical trials is rapidly increasing in order to quickly generate the data required for efficient patient healthcare. However, we have found that some of the pharmacology data published so far are somewhat disappointing, due to a lack of information permitting adequate comprehension of the dose-exposure and dose-effect relationships [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF][START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF] and the poor representativeness of data used for simulations of effective dosing regimens [START_REF] Garcia-Cremades | Optimizing hydroxychloroquine dosing for patients with COVID-19: An integrative modeling approach for effective drug repurposing[END_REF][START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF][START_REF] Yao | In Vitro Antiviral Activity and Projection of Optimized Dosing Design of Hydroxychloroquine for the Treatment of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)[END_REF].

For these reasons and on behalf of the Clinical Pharmacology Committee of the French agency for AIDS and viral hepatitis research (ANRS) and the Therapeutic Drug Monitoring and Treatment Personalization working group of the French Society of Pharmacology and Therapeutics (SFPT), we believe that there is an urgent need for clarifications and improvements in order to generate high-quality PK and PK-PD data for the drugs to be used for COVID-19 treatment.

Limitations of available PK and PK-PD data.

Few studies have already described PK and PK-PD in potential treatments for COVID-19 infection and only a small fraction of ongoing clinical trials have planned to do so. We were only able to find 18 out of the 1546 registered in clinicaltrials.gov (Table I) (www.clinicaltrials.gov, accessed on: May 18, 2020). The list presented here is representative of the current situation though not necessarily exhaustive. The situation is particularly striking for the commonly used repurposed drugs HCQ and LPV for which only 2 trials, one in children and another in adults, were designed to measure concentration data (Table I). PK analyses are preferentially planned for new drugs (i.e. monoclonal antibodies) not yet approved in another indication, such as to take advantage of COVID-19 infection to acquire PK data in humans. [START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF] and [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF] have reported HCQ blood or serum concentrations from COVID-19 patients. However, in these two papers, very little methodological, demographic, clinical, paraclinical and even dosing information that could help to better understand the dose-exposure relationship in COVID-19 patients is available. Perinel et al. have succinctly described a small population of intensive care unit (ICU) patients receiving 200 mg of HCQ three times daily without loading dose in which HCQ whole blood concentrations were measured [START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF]. A similar dosing regimen was used by [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF] but serum concentrations were determined instead. Apparently, for their study, they have reported the sum of HCQ and its metabolite concentrations. Whether the approximation used to quantify HCQ metabolite is valid [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF] and whether this metabolite is active on SARS-CoV-2 require further evaluation. Anyway, simply summing-up the concentrations of the active moiety and its metabolite cannot effectively contribute to comprehension of PK and PK-PD relationships since they may present different PK and/or PD properties and consequently misrepresent the adequate drug exposure [START_REF] Tett | Clinical pharmacokinetics of slow-acting antirheumatic drugs[END_REF]. These preliminary observations also raise the question of the selection of the adequate biological matrix (blood vs plasma or serum) for PK assessment of HCQ. Some of these papers [START_REF] Garcia-Cremades | Optimizing hydroxychloroquine dosing for patients with COVID-19: An integrative modeling approach for effective drug repurposing[END_REF][START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF][START_REF] Yao | In Vitro Antiviral Activity and Projection of Optimized Dosing Design of Hydroxychloroquine for the Treatment of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)[END_REF] have built dosing recommendations for HCQ based on PK models previously developed in other diseases (SEL, malaria) and clinical situations (steady-state, non-ICU patients), which raises questions about the relevance of such recommendations. Perinel et al. [START_REF] Perinel | Towards Optimization of Hydroxychloroquine Dosing in Intensive Care Unit COVID-19 Patients[END_REF] have graphically compared measured HCQ blood concentrations in COVID-19 patients to PK simulations obtained with a HCQ population PK model initially developed for rheumatoid arthritis (RA) patients [START_REF] Carmichael | Population pharmacokinetics of hydroxychloroquine in patients with rheumatoid arthritis[END_REF]. Results clearly demonstrate the inadequacy of a population PK model developed for stable chronic RA patients to describe actual PK data obtained in ICU patients. This point has also been confirmed by Martin-Blondel et al (Martin-Blondel et al., n.d.) in their attempt to describe HCQ plasma concentrations obtained in COVID-19 patients using a PK model initially developed for SLE patients. As a consequence, we believe that due to probably large PK differences affecting either clearance or volume of distribution as demonstrated with other anti-infective drugs, this type of PK model should not be used to simulate any dosing regimens for ICU patients [START_REF] Roberts | International Society of Anti-Infective Pharmacology and the Pharmacokinetics and Pharmacodynamics Study Group of the European Society of Clinical Microbiology and Infectious Diseases[END_REF]. Whether these models could be applied to less severe COVID patients, at the early stage of the disease, requires confirmation. It is likely that models developed in healthy volunteers may be appropriate in the prophylaxis setting. The population PK-PD model developed by [START_REF] Garcia-Cremades | Optimizing hydroxychloroquine dosing for patients with COVID-19: An integrative modeling approach for effective drug repurposing[END_REF] to describe viral decline and QTc prolongation after HCQ administration in COVID-19 patients, may theoretically enable achievement of more robust simulation results. For the PK part of the model, they used a population PK model initially developed using plasma concentrations obtained from healthy Korean volunteers and Korean patients receiving HCQ for treatment of vivax malaria [START_REF] Lim | Pharmacokinetics of hydroxychloroquine and its clinical implications in chemoprophylaxis against malaria caused by Plasmodium vivax[END_REF]. To demonstrate the adequacy of this model for COVID-19 patients, they used published serum concentrations obtained in non-severe and non-ICU COVID-19 patients drawn from the study by [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF] described above. Although some differences between serum and plasma concentrations cannot be ruled out [START_REF] Bergqvist | Distribution of chloroquine and its metabolite desethylchloroquine in human blood cells and its implication for the quantitative determination of these compounds in serum and plasma[END_REF], they showed that HCQ concentrations from COVID-19 patients fell within the lower range of expected population profiles, thereby suggesting that this model could be appropriate to describe the PK of HCQ, at least in non-ICU COVID-19 patients. However, one should note that the concentration data set used for this model validation expressed HCQ concentration as the sum of the active moiety and its metabolite [START_REF] Gautret | Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial[END_REF], consequently leading to systematic overestimation of actual exposure to HCQ. Other limitations of this paper apply to the PD parts of the model since viral kinetic data are from SARS-CoV-1 and the concentration/QTc prolongation data were obtained from a study with chloroquine. Limitations also apply to the study by [START_REF] Yao | In Vitro Antiviral Activity and Projection of Optimized Dosing Design of Hydroxychloroquine for the Treatment of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)[END_REF], who developed a PB-PK model to simulate different dosing regimens for HCQ in COVID-19 patients. Indeed, HCQ concentrations in lung were simulated by incorporating in their model blood-to-lung concentration data obtained by analyzing tissue homogenates after dosing in rats [START_REF] Mcchesney | Animal Toxicity and Pharmacokinetics of Hydroxychloroquine Sulfate[END_REF]. This approach may provide unreliable diffusion data in lungs [START_REF] Mouton | Tissue Concentrations: Do We Ever Learn? The Journal of antimicrobial chemotherapy[END_REF][START_REF] Nix | Antibiotic tissue penetration and its relevance: models of tissue penetration and their meaning[END_REF] and consequently unreliable simulations at the site of infection. Similar limitations of this lung model have been discussed by Fan et al [START_REF] Fan | Connecting hydroxychloroquine in vitro antiviral activity to in vivo concentration for prediction of antiviral effect: a critical step in treating COVID-19 patients[END_REF] and by Yeo et al [START_REF] Smith | Dosing will be a key success factor in repurposing antivirals for COVID-19[END_REF].

Similarly, a steady-state population PK model initially developed from HIV patients receiving standard dose LPV/r has been used to simulate LPV total and free concentrations for comparison with various virological endpoints, i.e. actual IC50 values for HIV-1, MERS-CoV and SARS-CoV [START_REF] Smith | Dosing will be a key success factor in repurposing antivirals for COVID-19[END_REF]. From a PK perspective, this model is probably unable to describe actual LPV concentrations currently measured in COVID-19 patients, which are as recently reported [START_REF] Gregoire | Lopinavir pharmacokinetics in COVID-19 patients[END_REF] much higher than those observed in HIV patients .

Furthermore, this recent paper showed that despite very high plasma LPV concentrations observed in patients, the free fraction, representing the active part of the drug, was not affected [START_REF] Gregoire | Lopinavir pharmacokinetics in COVID-19 patients[END_REF].

Concerning favipiravir, a recent comment proposes to evaluate higher doses for SARS-CoV-2 than in Ebola disease, based on both in vitro EC50 and PK simulations data associated with a close monitoring of adverse events and plasma concentrations [START_REF] Eloy | Dose rationale for favipiravir use in patients infected with SARS-CoV-2[END_REF]. However, favipiravir is a prodrug, and determination of tissue and intracellular exposure of the activate metabolite favipiravir-ribofuranosyl-5-triphosphate would be required in order to better characterize PK-PD relationship [START_REF] Du | Response to "Dose rationale for favipiravir use in patients infected with SARS-CoV-2[END_REF].

No such attempts have been made with remdesivir, yet. Furthermore, data on the diffusion of these drugs into the pulmonary tract are lacking.

Recommendations for improvements

We propose here to identify some of the major requirements that need to be addressed in designing PK and PK-PD studies in this era of COVID.

Pre-analytical and analytical requirements

Pre-analytical and analytical steps should be strictly controlled to guarantee the necessary accuracy of measured concentrations. For LPV, well-standardized sampling and assay procedures have been defined for the monitoring of HIV patients and we believe that these procedures may adequately apply to COVID-19 patients. By contrast, this is more confusing with HCQ since it could be analyzed either in whole blood [START_REF] Carmichael | Population pharmacokinetics of hydroxychloroquine in patients with rheumatoid arthritis[END_REF] or in plasma [START_REF] Morita | Population Pharmacokinetics of Hydroxychloroquine in Japanese Patients With Cutaneous or Systemic Lupus Erythematosus[END_REF][START_REF] Tett | Bioavailability of hydroxychloroquine tablets in healthy volunteers[END_REF]. In the area of anti-infective drugs, plasma may appear as the preferred matrix since free plasma concentration is in equilibrium with tissue concentrations and the putative site of action. Moreover, some PK-PD relationships have been commonly developed for antivirals using plasma concentration as a surrogate for the concentrations at the site of action [START_REF] Rizk | Pharmacokinetics and pharmacodynamics of once-daily versus twice-daily raltegravir in treatment-naive HIV-infected patients[END_REF]. However, as regards plasma, it appears that the pre-analytical step is of particular importance for chloroquine derivatives since early observations have shown that chloroquine could be released from blood cells, leading to overestimation of plasma concentrations [START_REF] Bergqvist | Distribution of chloroquine and its metabolite desethylchloroquine in human blood cells and its implication for the quantitative determination of these compounds in serum and plasma[END_REF]. Preliminary pre-analytical data suggest a similar pattern with HCQ when blood is not rapidly centrifuged and plasma separated. Moreover, HCQ blood-to-plasma ratio presents wide between-subject variability according to disease context and equilibrium distribution [START_REF] Morita | Population Pharmacokinetics of Hydroxychloroquine in Japanese Patients With Cutaneous or Systemic Lupus Erythematosus[END_REF][START_REF] Tett | A dose-ranging study of the pharmacokinetics of hydroxy-chloroquine following intravenous administration to healthy volunteers[END_REF] suggesting that extrapolation from one to another is not reliable. Pre-analytical treatment of plasma for remdesivir assay requires an acidic stabilization step (Gilead, personal communication). Regarding the analytical step, robust specific and sensitive assays targeting individual analytes (active moiety and metabolites) have already been developed and validated for clinical pharmacokinetics. In the event of matrix modification (blood vs. plasma), a partial validation could be considered (European Medicines Agency. Guideline on bioanalytical method validation, 2011). In addition, inter-laboratory comparisons should be organized, especially in case of multicenter PK studies.

Collection of covariates affecting dose-exposure relationships

In order to comprehensively understand and describe dose-exposure relationships, collection of accurate dosing information, demographic, clinical and para-clinical data is mandatory.

Regarding dosing information, galenic forms and routes of administration are important data.

It is common practice in the setting of an ICU, where the patient is unable to swallow solid oral dosage forms, to administer drugs through enteral feeding tubes. In this situation, the administration of crushed tablets in the feeding tube may have an impact on bioavailability, which may alter the dosing-exposure relationship of the studied drug. This has been demonstrated in the HIV setting where crushing tablets of LPV have led to a 45% decrease in AUC [START_REF] Best | Pharmacokinetics of lopinavir/ritonavir crushed versus whole tablets in children[END_REF]. Surprisingly, crushing tablet does not appear to have a significant impact on the oral bioavailability of HCQ (Sanofi, personal communication).

The dosing regimen may include a loading dose that is particularly important in view of more rapidly achieving effective concentrations. This point is critical for HCQ due to its very long half-life [START_REF] Lê | Rationale of a loading dose initiation for hydroxychloroquine treatment in COVID-19 1 infection in DisCoVeRy trial[END_REF][START_REF] Tett | Bioavailability of hydroxychloroquine tablets in healthy volunteers[END_REF]. A loading dose has also been suggested for LPV [START_REF] Smith | Dosing will be a key success factor in repurposing antivirals for COVID-19[END_REF] and is required for remdesivir (European Medicines Agency. Summary on remdesivir compassionate use. April 2020).

Knowledge of co-administered drugs is likewise of primary importance, particularly in the context of ICU, where many medications are involved or when complications such as secondary infections are present, requiring co-administration of antibacterials and antifungals [START_REF] Sanders | Pharmacologic Treatments for Coronavirus Disease 2019 (COVID-19): A Review[END_REF]. Co-medications are also frequent in patients particularly at risk for severe COVID-19, i.e. chronic disease patients such as obese and/or diabetic patients. Indeed, drug-drug interactions may occur at each step of the ADME process, leading to altered pharmacokinetics and, possibly, in increased variability of drug exposure. Divalent cations, such as calcium or magnesium, may interfere with the absorption of chloroquine [START_REF] Mcelnay | The effect of magnesium trisilicate and kaolin on the in vivo absorption of chloroquine[END_REF], thereby reducing its bioavailability. Similarly, LPV and HCQ, which undergo CYP mediated metabolism (Liverpool HIV group. Liverpool COVID-19 Interactions. https://www.covid19-druginteractions.org/, 2020), could be victims of drug-drug interactions affecting both bioavailability and clearance. Although remdesivir is also a substrate for CYP isoenzymes in vitro, its metabolism is likely to be predominantly mediated by hydrolase activity (Liverpool HIV group. Liverpool COVID-19 Interactions. https://www.covid19druginteractions.org/, 2020), but this does not stave off inter-nor intra-patient variability.

Anthropometric characteristics, which are currently involved in PK variability, should be collected as well. Body weight has been described as a significant covariate of HCQ clearance in Japanese patients suffering from systemic lupus erythematosus [START_REF] Morita | Population Pharmacokinetics of Hydroxychloroquine in Japanese Patients With Cutaneous or Systemic Lupus Erythematosus[END_REF].

Whether this relationship is valid in COVID-19 patients requires confirmation.

The use of specific organ support techniques [START_REF] Zhang | Response of Chinese Anesthesiologists to the COVID-19 Outbreak[END_REF] such as renal replacement therapy (hemodialysis, hemofiltration, hemodiafiltration) or ventilation support (high-flow oxygen therapy, noninvasive ventilation, invasive mechanical ventilation and extracorporeal membrane oxygenation, ECMO) should be reported in order to assess their impact on drug exposure. A significant impact of hemodialysis on the pharmacokinetics of HCQ [START_REF] Tett | Bioavailability of hydroxychloroquine tablets in healthy volunteers[END_REF] and LPV [START_REF] Gupta | The pharmacokinetics and pharmacogenomics of efavirenz and lopinavir/ritonavir in HIVinfected persons requiring hemodialysis[END_REF] is unlikely but requires confirmation since steady-state conditions are not achieved in COVID-19 patients. An impact of ECMO on chloroquine pharmacokinetics has been described [START_REF] Bagate | Early extracorporeal membrane oxygenation for cardiovascular failure in a patient with massive chloroquine poisoning[END_REF] but no convincing data exist yet for LPV [START_REF] Ghazi Suliman | The effect of veno-venous ECMO on the pharmacokinetics of Ritonavir, Darunavir, Tenofovir and Lamivudine[END_REF]. However, a significant impact of ECMO should be considered since several studies have shown altered PK profiles in this situation [START_REF] Ha | Evaluation of Altered Drug Pharmacokinetics in Critically Ill Adults Receiving Extracorporeal Membrane Oxygenation[END_REF]. These data need to be confirmed in COVID-19 patients, especially considering differences in steady-state attainment.

Importantly, severe COVID-19 is associated with a systemic hyper-inflammation state, the socalled cytokine storm, which is associated with highly elevated C-reactive protein (CRP) levels, hyperferritinaemia and increased cytokine levels (IL-1β, IL-2, IL-6, IL-17, IL-8, TNF and CCL2) [START_REF] Mcgonagle | The Role of Cytokines including Interleukin-6 in COVID-19 induced Pneumonia and Macrophage Activation Syndrome-Like Disease[END_REF]. The impact of inflammation on the pharmacokinetics of drugs has already been highlighted. For example, several studies have reported a positive association between CRP levels and voriconazole or tacrolimus concentrations [START_REF] Bonneville | Unexpected overdose blood concentration of tacrolimus: Keep in mind the role of inflammation[END_REF][START_REF] Gautier-Veyret | Inflammation is a potential risk factor of voriconazole overdose in hematological patients[END_REF]. The extent of increase in voriconazole concentration could be explained by downregulation of CYP isoenzymes by inflammatory stimuli leading to reduced metabolism [START_REF] Morgan | Impact of infectious and inflammatory disease on cytochrome P450-mediated drug metabolism and pharmacokinetics[END_REF]. Impact may also occur through altered expression of membrane transporters [START_REF] Seifert | Inflammation and pharmacokinetics: potential implications for HIV-infection[END_REF]. Preliminary data suggest a major decrease of LPV clearance in COVID-19 patients, highlighting the putative role of inflammation in these PK alterations.

Development of population PK and PK-PD models

Population approaches, which rely on PK-PD modeling, appear particularly appealing in this situation because these models can handle sparse data originating from different sources (ICU, non-ICU, dialysis patients) and different dosing regimens. Further, covariates as described above can be included in these population models to identify PK and pharmacodynamic (PD) variability factors. Finally, population-PK and PK-PD models have proven useful in the area of infectious diseases, providing a valuable tool to explore and predict efficient dosing regimens [START_REF] Jumah | Pharmacokinetic/Pharmacodynamic Determinants of Vancomycin Efficacy in Enterococcal Bacteremia[END_REF]. However, we believe that the development of population PK models specifically dedicated to acute COVID-19 patients is urgently needed due to major physiopathological differences with chronically ill HIV, RA and SEL patients. For that purpose, a collaborative multicenter study within clinical centers involved in clinical trials with these drugs may enable collection of the amount of data necessary to develop such a PK model. Such a collaborative project is ongoing within our groups.

Current dosing regimens proposed for these repurposed drugs are empirical since they have not been specifically developed for COVID-19 but rather for HIV, malaria and chronic inflammatory diseases. In order to develop more efficient dosing regimens, we also need reliable PD data to relate the concentration at the site of infection to viral sensitivity through the development of PK-PD models. Viral sensitivity could be determined through the measurement of EC50 or EC90 either in presence or in absence of human serum. Addition of human serum enables representation of the impact of drug protein binding on drug potency. EC50 is available for HCQ [START_REF] Liu | Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro[END_REF][START_REF] Yao | In Vitro Antiviral Activity and Projection of Optimized Dosing Design of Hydroxychloroquine for the Treatment of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)[END_REF], LPV [START_REF] Choy | Remdesivir, Lopinavir, Emetine, and Homoharringtonine Inhibit SARS-CoV-2 Replication in Vitro[END_REF], favipiravir [START_REF] Eloy | Dose rationale for favipiravir use in patients infected with SARS-CoV-2[END_REF], remdesivir [START_REF] Choy | Remdesivir, Lopinavir, Emetine, and Homoharringtonine Inhibit SARS-CoV-2 Replication in Vitro[END_REF] and other drugs such as ivermectin (Table II). But, as observed in Table II, translation into effective therapies may be challenging due to unfavorable in vivo PK properties. For example, ivermectin, total (bound and unbound) plasma concentrations are >45 times lower than the in vitro EC50. Since ivermectin is highly bound to plasma protein and its accumulation into human lung unknown [START_REF] Schmith | The Approved Dose of Ivermectin Alone Is Not the Ideal Dose for the Treatment of COVID-19[END_REF], the likelihood for ivermectin to reach IC50 at the site of action after the approved dose is low. The large variability reported in EC50 values between studies, mainly explained by differences in experimental conditions or mechanism of action, may also contribute to the difficulty of determining the accurate target concentration in vivo. For example, with HCQ, we observed high variability between EC50 values (i.e. for HCQ: from 0.72 μM ≈ 0.24 μg/mL to more than 15 μM ≈ 5 μg/mL) depending on MOIs, incubation duration, presence/absence of human serum, etc…, which may be of PK relevance for a drug presenting such a narrow therapeutic index. As reported in Figure 1, HCQ median plasma concentrations could be 1.5 to 15 times lower than the in vitro EC50 depending on the scheme of administration (and EC50 values. Moreover, the choice of experimental cells (Vero E6 cells or pulmonary epithelial cells) may be crucial for such respiratory viruses, although no consensus has yet arisen, and extrapolations from other viruses appear speculative (Smith et al., 2020). Finally, regarding the lack of specificity/sensitivity of the current virological marker (i.e viral load in the upper respiratory tract) to assess antiviral activity, the combination of multiple PD parameters could be of interest. The relationships between viral load and clinically relevant endpoints is mostly unknown. Indeed, many clinical endpoints could be considered: mortality, recovery, hospital discharge, respiratory failure, need for oxygenation or mechanical ventilation, time-to-intubation, hospitalization, duration of hospital-stay, ICU admission, time-to-improvement, severity scores (SOFA), etc. Selecting the best clinically relevant endpoint may depend upon several factors such as its clinical relevance, its measurement's reliability, disease severity and statistical considerations.

In order to assess drug penetration to the site of action, determination of drug concentration in bronchoalveolar lavage could be proposed as a surrogate for lung concentrations. Urea or albumin should be used as a marker of dilution to determine the volume of epithelial lining fluid recovered and samples should be preferably collected at multiple time points throughout the dosing interval [START_REF] Rodvold | Penetration of anti-infective agents into pulmonary epithelial lining fluid: focus on antibacterial agents[END_REF]. To our knowledge, corresponding data for COVID-19 treatment candidates are scarce. A case study in a single HIV patient [START_REF] Atzori | Detection of intrapulmonary concentration of lopinavir in an HIV-infected patient[END_REF] reported significant LPV concentration in BAL at steady-state. ELF and total plasma concentrations were 14.4 and 8.1 µg/mL, respectively. The corresponding ratio of ELF over total plasma concentrations was 1.8, suggesting an accumulation of LPV in ELF requiring confirmation in COVID-19 patients. Data of the same order of magnitude but with less accumulation were also reported by Boffito et al [START_REF] Boffito | Lopinavir measurement in pleural effusion in a human immunodeficiency virus type 1-infected patient with kaposi's sarcoma[END_REF].

Conclusion

Available PK and PK-PD studies suffer from severe limitations leading to unreliable conclusions, especially in term of dosing optimization. At this time, there is still no high-quality evidence to support the use of these repurposed drugs for the treatment of COVID-19.

Therefore, implementation of well-designed PK and PD studies specifically targeted to COVID-19 patients is urgently needed in order to increase our comprehension of the doseexposure-effect relationships of the repurposed drugs. Without these data, efficient evaluation and development of effective dosing regimens will remain difficult. Other scientific societies [START_REF] Baker | A call for the appropriate application of clinical pharmacological principles in the search for safe and efficacious COVID-19 (SARS-COV-2) treatments[END_REF] and organizations [START_REF] Hartman | Time to Step Up: A Call to Action For the Clinical and Quantitative Pharmacology Community to Accelerate Therapeutics for COVID-19[END_REF][START_REF] Rayner | Optimizing COVID-19 Candidate Therapeutics: Thinking Without Borders[END_REF] have also recently issued call to action for the appropriate application of clinical pharmacology principles in the search for COVID-19 treatments. Further, to accelerate the production of upto-date PK and PD data and the development of meaningful PK and PK-PD models, we also call for multi-institutional collaborative work and involvement of clinical pharmacologists in multidisciplinary research consortia. Lopinavir/ritonavir (400 lopinavir mg/100 mg ritonavir) will be administered every 12 h for 14 days in tablet form. For patients who are unable to take medications by mouth, the lopinavir/ritonavir (400 lopinavir mg/100 mg ritonavir) will be administered as a 5-ml suspension every 12 h for 14 days via a pre-existing or newly placed nasogastric tube.

Interferon ß1a will be administered subcutaneously at the dose of 44 µg for a total of 3 doses in 6 days (day 1, day 3, day 6).

Hydroxychloroquine will be administered orally as a loading dose of 400 mg twice daily for one day followed by 400 mg once daily for 9 days. The loading dose of hydroxychloroquine through a nasogastric tube will be increased to 600 mg twice a day for one day, followed by a maintenance dose of 400 mg once a day for 9 days

Figure 1 .

 1 Figure 1. Observed hydroxychloroquine plasma concentrations (median) after 3-5 days of treatment in patients treated for SARS-CoV-2 infection depending on the scheme of administration: 400mg QD after a loading dose of 400mgx2 on day 1 and 200mg TID without any loading dose (C. Solas, internal data from the Pharmacokinetics and Toxicology Laboratory). Solid line represent in vitro EC50 described against SARS-CoV-2 at T48h post incubation.

Table I :

 I Pharmacological issues relative to COVID-19 infection in clinical trials : 18 studies found for: Pharmacokinetic | COVID-19 579 580

		Intravenous Infusion to Subjects With	events respiratory failure [Time Frame:	Hydrochloride, Spironolactone,
		Yellow Fever or COVID-19	Day 14]	Terbutaline, Tranexamic acid
	Clinical trial Identifier NCT04345614 NCT04350736 NCT03648372	Official Title A Randomized Controlled Open-Label Study of CM4620 Injectable Emulsion in Patients with Severe COVID-19 Pneumonia A Phase 1, Double-blind, Randomized, An Open Label, Dose-Escalation, Phase I Study to Evaluate the Safety, Tolerability and Pharmacokinetics of TAK-981 in Adult Patients With Advanced or Metastatic Solid Tumors or Relapsed/Refractory Hematologic Malignancies and in a Subset Placebo-controlled, Sponsor-open, SAD and MAD Study in Healthy Subjects to Evaluate Associated With COVID-19 TD-0903, a Potential Treatment for ALI the Safety, Tolerability, and PK of Inhaled	Primary outcome Safety, efficacy, and the PK profile of CM4620-IE in patients with severe COVID-19 pneumonia. 2-change in laboratory Safety and Tolerability of MAD of parameters COVID-19 proof of concept: Once the Safety Lead-in is complete and a TAK-981 dose and regimen is selected by the Safety Monitoring Committee (SMC), the randomized COVID-events treatment emergent adverse Number and severity of TD-0903: Adverse Events [ Time Frame: Day 1 to Day 14 ]	Experimental drug CM4620-Injectable Emulsion Voriconazole, Zolpidem TD-0903 TAK-981 Intravenous infusion Dextroamphetamine /Amphetamine -Extended Release, Azithromycin Chloroquine, Hydroxychloroquine, Lopinavir/Ritonavir, Ribavirin, Tocilizumab inhalation Study drug to be administered by	Pharmacological M4620-IE serum concentration Plasma PK parameters of SAD TAK-981 COVID-19 safety lead-in: TAK-981, intravenously, administered as 60 minute-infusion, once on Days 1 and 4. The starting dose of TAK-981 will be 60 milligram (mg). and MAD TD-0903
	NCT04357613 NCT04346628 NCT04371640 NCT04392219 NCT04369469 NCT04379271	A RANDOMIZED NON-COMPARATIVE PHASE 2 PILOT STUDY TESTING THE VALUE OF IMATINIB MESYLATE AS AN EARLY TREATMENT OF COVID-19 DISEASE IN AGED HOSPITALIZED PATIENTS. A Phase 2 Randomized, Open Label Study of Oral Favipiravir Compared to Standard Supportive Care in Subjects With Mild COVID-19 With Coronavirus Disease 2019 A Randomized, Double-Blinded, Placebo-Controlled Trial Evaluating the Virological Efficacy, Safety, Tolerability, Pharmacokinetics, and Pharmacodynamics of Sirolimus Adjuvant Therapy in Patients With COVID-19 A Randomized, Double-Blind, Placebo-Controlled, First-in-Human Study Designed to Evaluate the Safety, Tolerability, and Pharmacokinetics of EIDD-2801 Following Oral Administration to Healthy Volunteers A Phase 3 Open-label, Randomized, Controlled Study to Evaluate the Efficacy and Safety of Intravenously Administered Ravulizumab Compared With Best Supportive Care in Patients With COVID-19 Severe Pneumonia, Acute Lung Injury, or Acute Respiratory Distress Syndrome A Prospective, Multi-Center, Randomized, Placebo-Controlled, Double-Blinded Study to Evaluate the Efficacy, Safety and Tolerability of IMU-838 as Addition to	To evaluate the benefit of early imatinib therapy to prevent severe COVID-19 disease in hospitalized aged patients. [Time Frame: 30 days] Time until cessation of oral shedding of SARS-CoV-2 virus [ Time Frame: Up to 28 days ] 19 proof of concept will begin with participants randomized to Arm A: COVID-19 standard of care (SOC), or Arm B COVID-19 SOC + TAK-981. Change in SARS-COV-2 viral burden from baseline to day 7 of treatment [Time Frame: Baseline, and days 1, 2, 3, 4, 5, 6, & 7 post-dose for all patients] Safety and Tolerability of Single Ascending Dose (SAD) of EIDD-2801 (Part 1) and Multiple Ascending Dose (MAD) of EIDD-2801 (Part 3): Adverse Events [ Time Frame: From screening Survival (based on all-cause mortality) at Day 29 [ Time Frame: Baseline, Day 29 ] through study completion, up to 15 days ] Number and severity of treatment emergent adverse events Proportion of patients without any need for INV until end-of-study (EoS) [ Time Frame: Throughout the Study (Day 0 to	Imatinib 800mg/d during 14days Favipiravir administered orally, 1800 mg on the first dose (day 1) followed by 800 mg twice daily for the next 9 days (days 2-10). Sirolimus + standard medical care Day 1: 10mg Days 2-7: 5mg Intervention: Drug: Sirolimus 1 MG/ML Part 1: Subjects will be randomized to receive a single oral dose of EIDD-2801 or Placebo. Part 3: Subjects will be randomized Weight-based doses of ravulizumab will be administered intravenously on Days 1, 5, 10, and 15. Other Names: Ultomiris to receive twice daily dosing either ALXN1210 EIDD-2801 or Placebo.	To evaluate plasma levels of imatinib [Time Frame: 14 days] Imatinib trough level Cmax of favipiravir [Time Frame: Days 1 and 10 (samples taken 30 minutes prior to and 1 hour following favipiravir administration)] PK of sirolimus PK parameters of ravulizumab Pharmacokinetics (PK) of EIDD-2801 when given as Single Doses (Part 2): Maximum observed concentration Cmax [ Time Frame: Day 1 through Day 18 ] Multiple pharmacokinetic variables of EIDD-2801 will be assessed and may include, but are not limited to: Maximum observed concentration Cmax
	NCT04358549 NCT03891420 NCT04320615 NCT04363736 NCT04380688 NCT04346199	Open Label, Randomized, Controlled Phase 2 Proof-of-Concept Study of the Use of Favipiravir v. Standard of Care in Hospitalized Subjects With COVID-19 A Phase 1b Double-blind, Placebo-controlled, Dose-ranging Study to Evaluate the Safety, Pharmacokinetics , and Anti-viral Effects of Galidesivir Administered Via A Randomized, Double-Blind, Placebo-With COVID-19 A Phase-II, Open-Label, Randomized, Multicenter Study to Investigate the Pharmacodynamics, Pharmacokinetics, Safety, and Efficacy of 8 mg/kg or 4mg/kg Supportive Care in Subjects Hospitalized Moderate to Severe COVID-19 Pneumonia With Best Supportive Care Versus Best Intravenous Tocilizumab in Patients With of the Efficacy and Safety of Acalabrutinib A Phase 2, Open Label, Randomized Study With COVID-19 Supportive Care in Subjects Hospitalized A Phase 2, Open Label, Randomized Study of the Efficacy and Safety of Acalabrutinib With Best Supportive Care Versus Best Disease 19 Therapy, in Patients With Coronavirus Investigator's Choice of Standard of Care	Time to viral clearance [Time Frame: Day 29]; To determine the effect of favipiravir + SOC v. SOC on viral clearance of COVID-19 as measured by nasopharyngeal and oropharyngeal sampling Number of subjects with 1-treatment emergent adverse events and serious adverse Clinical Status Assessed Using a Subject alive and free of dose] [Time Frame: 28 days after last Concentration of C-Reactive 7] and Serious Adverse Events Occurrence of Adverse Events Protein (CRP) [Time Frame: Day SARS-COV-2 viral burden will be quantified for both arms using a qRT-PCR Subject alive and free of respiratory failure [ Time Frame: Day 14 ] Day 28) ]	Day 1: favipiravir 1800 mg BID plus Standard of Care (SOC) Days 2-14: 1000 mg BID plus SOC. For subjects Days 2-14: 800 mg BID plus SOC Galidesivir IV infusion Participants will receive 1 mg/kg vs 4 mg/kg in addition to standard-of-care treatment. Acalabrutinib-administered orally Participants will receive intravenous (IV) tocilizumab (TCZ) at a dose of 8 with Child-Pugh A liver impairment: Acalabrutinib-administered orally	Exposure of galidesivir as measured by plasma Serum Concentration of TCZ dose] concentrations [Time Frame: 28 days after last metabolite ACP-5862 acalabrutinib and its active Plasma PK parameters of marker of TCZ PK Protein (CRP) as surrogate last dose ] Concentration of C-Reactive [ Time Frame: 28 days after Plasma PK of favipiravir metabolite ACP-5862 acalabrutinib and its active Plasma PK parameters of

ACKNOWLEDGMENTS § ANRS-AC43 Clinical Pharmacology and SFPT Therapeutic Drug Monitoring and Treatment Personalization Committees: Zoubir Djerada, Matthieu Grégoire, Florian Lemaitre, Damien Montange, Patrice Muret, Jean-Marc Treluyer

We wish to thank Jeffrey Arsham, an American medical translator, for editing the English of our original manuscript.

Declarations of interest for NV, SB, CS, MCG, RG, RG and CS: None GP has received travel grants, consultancy fees, honoraria, or study grants from various pharmaceutical companies, including Gilead Sciences, Janssen, Merck and ViiV Healthcare. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

 , 2, 3, 4, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, and 49, andThe prescribing of drugs to children is not part of this protocol. Participants will receive DOIs as prescribed by their treating provider. Hydroxychloroquine EC50 = 0.72 µM (0.24 µg/mL) (Yao et al., 2020) EC50 = 4.51 µM (1.47 µg/mL) [START_REF] Liu | Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro[END_REF] EC50 = 4.17 µM (1.35 µg/mL) [START_REF] Touret | In vitro screening of a FDA approved chemical library reveals potential inhibitors of SARS-CoV-2 replication[END_REF] 200mg TID MD without LD: Ctrough = 0.09 +/-0.01 (H8 on day 2) to 0.19 +/-0.06 (H8 on day 6) µg/mL 400mgx2 LD on day 1 followed by 400x1 MD: Ctrough = 0.09 +/-0.12 (H12 on day 2) to 0.13 +/-0.14 (H12 on day 6) µg/mL (Martin-Blondel et al., n.d.) Favipiravir EC50 = 9.4 µg/mL EC50 = 40-80 µg/mL [START_REF] Eloy | Dose rationale for favipiravir use in patients infected with SARS-CoV-2[END_REF] EC50 > 100 µM (>15.71 µg/mL) [START_REF] Choy | Remdesivir, Lopinavir, Emetine, and Homoharringtonine Inhibit SARS-CoV-2 Replication in Vitro[END_REF] EC50 > 500 µM (>78.55 µg/mL) [START_REF] Jeon | Identification of antiviral drug candidates against SARS-CoV-2 from FDA-approved drugs[END_REF] Ebola patients 6000 mg on day 0 followed by 1200 mg BID for 9 days : Median (range) Ctrough = 46.1 (2.3 -106.9) µg/mL on day 2 and 25.9 (0 -173.2) µg/mL on day 4 [START_REF] Nguyen | Favipiravir pharmacokinetics in Ebola-Infected patients of the JIKI trial reveals concentrations lower than targeted[END_REF] Ivermectin EC50 = 2 µM (1750 ng/mL) [START_REF] Caly | The FDA-approved Drug Ivermectin Inhibits the Replication of SARS-CoV-2 in Vitro[END_REF]

Onchocerciacis patients

Single dose of 150 µg/kg PO: Cmax = 38.2 +/-5.8 ng/mL [START_REF] Okonkwo | Protein binding and ivermectin estimations in patients with onchocerciasis[END_REF] *Remdesivir is a prodrug rapidly converted to a circulating monophosphate nucleoside analog


Julien Joseph
email: julien.joseph@ens-lyon.fr

Fernanda Santibáñez
email: fsantibanez@unrn.edu.ar

María Fabiana Laguna
email: lagunaf@cab.cnea.gov.ar

Guillermo Abramson
email: abramson@cab.cnea.gov.ar

Marcelo N Kuperman
email: kuperman@cab.cnea.gov.ar

Lucas A Garibaldi
email: lgaribaldi@unrn.edu.ar

María Fabiana Laguna

A spatially extended model to assess the role of landscape structure on the pollination service of Apis mellifera

Keywords: pollination, Apis mellifera, competition, landscape structure, ecological intensification

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

Introduction

The intensive agricultural production systems developed during the twentieth century are facing multiple crises, mostly due to the degradation of the natural ecosystems they promote [START_REF] Newbold | Global effects of land use on local terrestrial biodiversity[END_REF][START_REF] Wilting | Quantifying biodiversity losses due to human consumption: a global-scale footprint analysis[END_REF][START_REF] Grab | Agriculturally dominated landscapes reduce bee phylogenetic diversity and pollination services[END_REF]. One of these crises is the fast decline of the populations of pollinators, and the consequent decrease in the quality of the pollination service [START_REF] Potts | Global pollinator declines: trends, impacts and drivers[END_REF][START_REF] Vanbergen | Threats to an ecosystem service: pressures on pollinators[END_REF]. The reduction of floral diversity requested by intensive monoculture particularly endangers wild pollinators, but also weakens managed species [START_REF] Potts | Global pollinator declines: trends, impacts and drivers[END_REF][START_REF] Vanbergen | Threats to an ecosystem service: pressures on pollinators[END_REF][START_REF] Potts | Safeguarding pollinators and their values to human well-being[END_REF][START_REF] Klein | Importance of pollinators in changing landscapes for world crops[END_REF]. Pollinators are not only essential for the functioning of terrestrial ecosystems, but they are also fundamental to our survival [START_REF] Potts | Safeguarding pollinators and their values to human well-being[END_REF][START_REF] Garibaldi | Policies for ecological intensification of crop production[END_REF]. Even if recent studies demonstrate the importance of wild bees, the managed bee Apis mellifera holds a particularly important place in the pollination service [START_REF] Garibaldi | Wild pollinators enhance fruit set of crops regardless of honey-bee abundance[END_REF][START_REF] Hung | The worldwide importance of honey bees as pollinators in natural habitats[END_REF].

In this context, countermeasures are constantly being developed to reach a better agriculture: productive, sustainable and healthy for the consumer [START_REF] Garibaldi | Policies for ecological intensification of crop production[END_REF][START_REF]UNEP Decisions Adopted by the Conference of the Parties to the Convention on Biological Diversity at its Fifth Meeting (UNEP/CBD/COP/5/23/Annex III), Decision V/5[END_REF][START_REF] Foley | Solutions for a cultivated planet[END_REF]. These include the structuring of landscapes to improve the ecological services provided by the biological actors of the crops' growth and reproduction, such as pollinators [START_REF] Garibaldi | Policies for ecological intensification of crop production[END_REF]. The experimental approach to study the effect of landscape structure on pollination quickly reaches limits in terms of the quantity of available experimental farms and of reproducibility of experiments. An approach in terms of mathematical modeling is particularly fit in such cases, since it allows the fast exploration of a variety of scenarios.

Landscape configurations found in agricultural exploitations around the world have arisen mostly from the experience of the growers, and are usually designed to increase the size of the cultivated area [START_REF] Macdonald | Farm size and the organization of U.S. crop farming[END_REF]. The treatments and management tested on the field or in silico to improve the pollination service usually consist in the addition of seminatural habitat-such as flower strips or grassy bands-mostly on the margins of the crops [START_REF] Hussler | Pollinator population size and pollination ecosystem service responses to enhancing floral and nesting resources[END_REF][START_REF] Scheper | Local and landscape-level floral resources explain effects of wildflower strips on wild bees across four European countries[END_REF][START_REF] Jnsson | Sown flower strips in southern Sweden increase abundances of wild bees and hoverflies in the wider landscape[END_REF]. But little is known about the global effects of the topology of the spatial arrangements of habitat in agricultural systems [START_REF] Martin | The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe[END_REF].

Many models have been developed in order to understand the needs and implications of sustainable pollination. Some of them consider a homogeneous environment and two interacting populations (plant and pollinator) [START_REF] Montoya | Trade-offs in the provisioning and stability of ecosystem services in agroecosystems[END_REF][START_REF] Cropp | Simulating eco-evolutionary processes in an obligate pollination model with a genetic algorithm[END_REF][START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF]. They generally do not include the spatial dimension, ignoring the effect of plants distribution in space or the distance between plants and pollinators' nests. Nevertheless, such models have provided an understanding of many characteristics of the plant-pollinator interactions in terms of the impact of the presence of semi-natural habitat [START_REF] Montoya | Trade-offs in the provisioning and stability of ecosystem services in agroecosystems[END_REF], or the coevolutive dynamics between plant and pollinator [START_REF] Cropp | Simulating eco-evolutionary processes in an obligate pollination model with a genetic algorithm[END_REF]. Besides, [START_REF] Croft | Evaluating European Food Safety Authority protection goals for honeybees ( Apis mellifera): What do they mean for pollination? Integrated environmental assessment and management[END_REF] [START_REF] Croft | Evaluating European Food Safety Authority protection goals for honeybees ( Apis mellifera): What do they mean for pollination? Integrated environmental assessment and management[END_REF], motivated by the growing concern on the abrupt and unexplained collapse of honeybee colonies, presented a model to evaluate the impact of pesticides and the regulation of their use on the sustainability of a colony. The authors concluded that the decrease of the colony size is a much better indicator of reduced ecoservice than an increased mortality of brood or foragers. While the work is focused on a very detailed description of the population dynamics of the individuals within the colony, with special attention on the effect of pesticides in three development stages, it does not incorporate details about the interaction between the bees and the resource. As a general feature, these models are not suitable to study the effect of landscape configurations, as the lack of a spatial dimension conceals the importance of a patterned distribution of resources, where boundaries and heterogeneous distributions can play extremely relevant roles, as will be discussed in the present work.

It is fair to mention, regarding the inclusion of a spatial dimension, several spatially explicit models that also incorporate seasonal aspects of plants and pollinators dynamics. They are able to predict the probable distribution of the pollinators on a large scale according to the resource distribution [START_REF] Lonsdorf | Modelling pollination services across agricultural landscapes[END_REF][START_REF] Olsson | Modeling pollinating bee visitation rates in heterogeneous landscapes from foraging theory[END_REF]. They can also infer the impact of some forms of management or treatments of a particular landscape on pollination [START_REF] Hussler | Pollinator population size and pollination ecosystem service responses to enhancing floral and nesting resources[END_REF].

In particular, the BEEHAVE family of models has addressed some questions on pollination using a spatially explicit approach [START_REF] Becher | BEEHAVE: a systems model of honeybee colony dynamics and foraging to explore multifactorial causes of colony failure[END_REF][START_REF] Becher | BEESCOUT: A model of bee scouting behaviour and a software tool for characterizing nectar/pollen landscapes for BEEHAVE[END_REF][START_REF] Becher | Bumble-BEEHAVE: a systems model for exploring multifactorial causes of bumblebee decline at individual, colony, population and community level[END_REF]. These models addressed the impact of a great number of stressors on pollinators at different scales (individuals, colony, population), while the model we present here focuses on the effect of landscape structure on flower visit rates. Bumble-BEEHAVE is the first model to include competition within and between hives and has shown its capacity to estimate at large scale the flower visit rates [START_REF] Becher | Bumble-BEEHAVE: a systems model for exploring multifactorial causes of bumblebee decline at individual, colony, population and community level[END_REF]. The study of this last aspect is the nucleus of another very interesting model, which is the one underneath the software tool BEESCOUT [START_REF] Becher | BEESCOUT: A model of bee scouting behaviour and a software tool for characterizing nectar/pollen landscapes for BEEHAVE[END_REF]. The model implemented in that tool points to examine the exploratory patterns of bees considering the particular structure of the landscape. It includes details about the detection capabilities of bees and propose several search strategies. The aim of this work is to contribute to the sim-ulation of realistic scenarios of colony growth and death in response to environmental and resource changes. In particular, the implementation of the foraging site election in these models leads to a systematic overchoosing of the most resourceful sites inside of a chosen "metapatch" of the grid. The sites with few resources are also excluded from the options of the bumblebees. Even if these approximations are useful to save computation time and are reasonable at some level to compare quite different landscapes, suboptimal choices due to partial knowledge of the surroundings and to communication failure between scouts and foragers have been shown to increase the dispersion of the visit and improve the pollination service [START_REF] Okada | Error in the honeybee waggle dance improves foraging flexibility[END_REF]. Our model uses parallel computation on GPUs (Graphics Processing Units) to avoid these approximations, and propose an adjustable site election modality allowing for mistakes. The scale at which our model operates is also finer (10 m against 25 m for Bumble-BEEHAVE) allowing to explore the question of the size of patches on the pollination efficiency. A comparison of our model to previous foraging models is shown in Table 1.

All the models cited above only consider a resource dynamic that is not affected by pollination [START_REF] Hussler | Pollinator population size and pollination ecosystem service responses to enhancing floral and nesting resources[END_REF][START_REF] Lonsdorf | Modelling pollination services across agricultural landscapes[END_REF][START_REF] Olsson | Modeling pollinating bee visitation rates in heterogeneous landscapes from foraging theory[END_REF][START_REF] Becher | BEEHAVE: a systems model of honeybee colony dynamics and foraging to explore multifactorial causes of colony failure[END_REF][START_REF] Becher | BEESCOUT: A model of bee scouting behaviour and a software tool for characterizing nectar/pollen landscapes for BEEHAVE[END_REF][START_REF] Becher | Bumble-BEEHAVE: a systems model for exploring multifactorial causes of bumblebee decline at individual, colony, population and community level[END_REF]. This is the second feature present exclusively in our model: a feedback of the resource in response to pollination. Even if, in general, nectar recuperation happens within a day [START_REF] Krlevska | Research on nectar-bearing capacity of apples[END_REF], a visited patch is likely to offer less resource on days following a visit because successfully pollinated flowers will stop producing pollen and nectar, and also because some of the flowers can be damaged by over-visitation [START_REF] Sáez | Extremely frequent bee visits increase pollen deposition but reduce drupelet set in raspberry[END_REF]. This has the effect of increasing the competition between hives due to a local short-term resource depletion.

Regarding the colony dynamics, our model uses the previously established growth model for Apis mellifera of [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF] [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF], combined with a novel foraging model which benefits from recent advances in parallel computation in order to infer the spatial distribution of visits to flowers from a fine resolution map of land cover compatible with a Geographic Information System (GIS). We consider only managed hives of A. mellifera whose wellbeing only depends on the available resource of their surroundings. We chose to model only managed hives because bee-keeping is a worldwide spread practice and the results we can produce will be applicable to more landscapes. In order to isolate the effects of the landscape structure we also compare the results obtained with the spatially explicit model to a similar mean-field model that assumes a uniform distribution of the different types of land covers.

To summarize, in this work we propose a step towards a theoretical foundation for agricultural landscape design, studying the effects of general compositional and structural features, such as the fraction of natural habitat, size of patches, and meshing of the patches of different crops on artificial maps, in order to be able to design new fields that optimize the pollination by A. mellifera. We intend to assess a number of relevant questions, such as: What is the influence of different configurations of natural habitats and crops (considering contrasting edge densities) on the dynamics of the honey-bee population, honey production, and crop flower visitation? How does the number of hives (as a management tool) interact with the effects of the configuration of natural habitats and crops through local competition? And also, what is the influence of different crop blooming periods? In the next section we describe the conceptual and mathematical details of our model, followed by results, a specific mean-field model to contrast them, and a final discussion.

Spatially explicit model

In this section we present the spatially explicit model. It is followed in the next section by a description of the mean field model, to be able to discuss the agreements and differences between them later on.

Landscape maps and resource dynamics

The maps used by the model are grids readable by GIS software, composed of a discrete number of habitat categories. Each habitat has a resource carrying capacity, a daily resource renewal rate and a blooming period. In this study we will consider that the carrying capacity is constant during the blooming period, and equal to zero the rest of the year, but a temporal distribution of the value of the carrying capacity can also be provided. The resource consists of both nectar and pollen considered together in the present model. Time evolves in discrete steps, with a time unit of one day.

When the resource is depleted a fraction of it disappears, proportionally to the number of visits, and it is renewed during the following days. Let K be the carrying capacity and R t the value of the resource before a number V t of visits on day t. Then, the resource of a patch evolves as:

R t+1 = R t -pV t + r(K -(R t -pV t )), (1)
where p is the fraction of flowers which will cease to produce nectar or pollen after the visits of bees (i.e. the success of pollination plus the damaged flowers), and r is a constant renewal rate. This equation reflects the decrease of the resource due to the daily visits and the recovery dynamics limited by saturation.

For the present study we created random maps of 3 km by 3 km, with each cell measuring 10 m by 10 m. This scale is fine enough to draw realistic maps relevant to the size of A. mellifera foraging sites (around 30 m by 30 m, i.e. 9 cells), and it is large enough to be interesting regarding their flight range, such that not all hives can forage on the entire map. These maps are composed of two habitats: a natural one whose blooming period is all spring and summer due to the great diversity of flowers it contains, and a monoculture with a 35 days blooming period. We used maps with different fractions of natural habitat and of edge density (number of cells of a given habitat which are at the border, divided by the total number of cells of the habitat). To build them, we used an algorithm of nucleation in 2D [START_REF] Kashchiev | Nucleation: Basic theory with applications[END_REF], varying the number of patches of natural habitat and the number of nuclei from which they grow. More nuclei mean smaller and more dispersed patches, and thus a higher edge density. The value of the resource is the same in the natural habitat and in the cultivated crop, so that the bees do not have any preference when the flowering periods overlap. This is a particular choice, used as a first step in this analysis, but can be easily relaxed when studying specific examples of crops and natural habitats.

The hives

Each hive is characterized by its position on the map, its population of bees and its quantity of honey. When the population of a hive falls below a threshold (one thousand bees), the hive is considered dead and removed from the list of hives.

Given that our main interest is the estimation of the number of visits to the crops, the only seasons considered in this study are spring and summer, when the crops bloom. Managed hives are usually well cared for during autumn and winter [START_REF] Furgala | Fall management and the wintering of productive colonies[END_REF] and we did not model their dynamic during this season. We only considered that at the beginning of a new foraging season, approximately one third of the bees of each hive has survived.

The dynamic of the hive is based on the one described by [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF] [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF]. The only difference is that, while in that model the food increases proportionally to the number of foragers, in our model it increases with the daily gain of the foragers, G t . The new equation for the food dynamic is the following:

F t+1 = F t -γ A (n f + n w ) -γ B n b + G t , (2)
with F t the quantity of food on day t, γ A the consumption rate of workers and foragers, γ B the consumption rate of broods (as in [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF]), and n f,w,b are the number of foragers, workers and broods, respectively. Apis mellifera is faithful to their harvesting sites: when they choose a site, they forage on it all day sending a great number of bees, sometimes until the resource of nectar is completely drained [START_REF] Rollin | Impacts of honey bee density on crop yield: A meta-analysis[END_REF]. In the model, each day a hive chooses a number of foraging sites with a size of 30 m by 30 m (9 cells), where it will send a fraction of foragers to harvest. The daily gain of the hive is the sum of the gains that each "squad" of foragers could obtain from the foraged site:

G t = M (Bt) i=1 g i,t (x, y), ( 3
)
where M is the number of squads which can forage from a hive of size B t (i.e. an integral fraction of the number of bees B t ), and g i,t (x, y) is the local gain of squad i at the foraging site centered at (x, y).

Apis mellifera choose their harvesting sites at the beginning of each day relying on relevant information provided by scouts [START_REF] Von Frisch | The dance language and orientation of bees[END_REF]: the resource at the sites and their distance from the nest. In the model, we implement this by choosing harvesting sites at random, at the beginning of each day, according to a specific distribution probability as follows. The probability to choose site (x, y) is defined as depending on R t (x, y), which is the total resource in the 9-cells Moore neighborhood centered at (x, y), and the flight cost to reach them, f (d):

P (x, y) =    (Rt(x,y)f (d)) γ (u,v)∈r f (Rt(u,v)f (d)) γ , if d < r f , 0, otherwise, (4)
where d is the Euclidean distance from the hive to (x, y) and r f is the range of flight of the bees. The cost f (d), normalized to the interval [0, 1], is an affine decreasing function of d, with f (0) = 1 and f (r f ) = 0. The exponent γ represents the knowledge that the hive has of the resource on its surroundings; if γ = 0, the hive chooses its foraging sites uniformly at random, and when γ → ∞, the hive systematically chooses the site with the highest harvestable resource. The precise value of this parameter is less important than the general shape of the function, so there is some freedom to choose it within the boundaries of biological significance. In our studies we used γ = 3, which gives a realistic behavior for the scouts, based on the shape of the preference function.

Once the harvesting sites are chosen, the hive sends a fraction ν i of foragers to each site i, also depending on resource and distance. This parameter ν i is a weight factor so that the more harvestable resource there is on a site, more bees are sent to it:

ν i = R t (x, y)f (d) chosen (x,y) R t (x, y)f (d) . ( 5
)
As a result of this election process, several hives can choose the same sites and thus compete for the resource. Note that the probability for hives to choose the same site is higher when the resource is globally low on their surroundings, and thus the competition is stronger.

We implement this possibility of hive competition by sorting each day a random harvesting order τ for each hive. Sorting this order every time step averages the effect of the chosen order. Let R τ t (x, y) correspond to the value of the resource after the turn of the hive which feeds at round τ . The actual harvestable resource for the next hive at the site is the product of R τ t (x, y) and f (d). Let c correspond to the maximum that a forager can carry back if the distance to the nest is null. If the resource is in excess, the bees sent to the site harvest the maximum value they can carry during the day, cf (d) ν k n f , but if the resource is lower than this value, they take it all. The daily local gain of hive i at site (x, y) can then be described as follows:

g i,t (x, y) = max (cf (d) ν i n f , R τ t (x, y)) . (6)
In the results we will also discuss the total number of visits, which is a proxy for the pollination efficiency. The number of visits to a site on one day is the sum of the visits made by all the visiting hives. In turn, the visits made by each of these hives is the product of the number of bees that flew to the site (N trips ), the average number of flowers visited in one trip (N f lower,trip ) and the average number of trips made by one bee:

V t = chosen sites N trips N f lower,trip ν i n f . (7)
Appropriate values are used for these parameters to compute the pollination service, as will be discussed below.

Setting the parameters

Since our main interest is the impact of landscape structure on pollination, all the parameters which concern the hives dynamic have been kept constant in the simulations (both in the spatially explicit model and in the mean field presented below). All the parameters used for the colony model have been kept from the study by [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF] [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF].

The parameters of the foraging model have been estimated from real data, such as the number of chosen sites and the number of trips per bee per day, and have been set within the range of observed values. A summary of the parameters and their values is provided in Table 2.

The average number of trips per day and the number of flowers visited have been chosen as the overall daily number of trips of a hive without distinguishing pollen trip from nectar trip.

The maximum distance of foraging has been set to 3 km so that bees can forage on almost all the map. It is known that bees can make trips farther than 3 km [START_REF] Abou-Shaara | The foraging behaviour of honey bees, Apis mellifera: a review[END_REF] but these trips cannot be seen on the map as currently implemented. The average foraging distance of the bees is not a parameter of the model, but a result of the foraging process, and varies with the exponent γ of the preference function and with the distribution of the resource.

The number of visits necessary to pollinate a flower is a parameter difficult to come by as it is highly dependent on the type of flower, and also depends on the scale considered in a nontrivial way. It is not available in the literature for the scale considered here, so we had to adjust it to observe a reasonable behavior of the resource depletion.

A second parameter which is not available in the literature is the size of a foraging squad. Nevertheless, as shown in the sensitivity analysis in the Supplementary Material, its value does not influence much the outputs of the model.

All the codes are written in Python using the CUDA module for GPU computation.

Mean-field model

Let us now describe a mean-field model of the system, which assumes a uniform distribution of the different types of land covers.

The dynamic of the hives during autumn and winter is the same as in the spatially explicit model: bees stop their foraging activities to enter in dormancy. For spring and summer, the mean-field equations for the number of bees in each hive and for the quantity of honey stored in the hive are also the same as in the spatially explicit model (from [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF]). However, the equation for the daily gain G t differs from Eq. ( 3), as now it is proportional to the number of foragers and to the available resource within the flight range of the hive, but has no dependence on the spatial distribution of the resource. In the following we discuss its dynamic in more detail.

As in the spatially explicit model, we consider that the hives are surrounded by two types of habitats: one is a monoculture relying on pollination, while the other is a natural habitat with greater flower diversity, such as a forest or even a semi-natural habitat. The monoculture has a determined blooming period, whereas the natural habitat is flowering during all spring and summer. We assume that during the blooming period the resource is in excess, and call R max the maximum that a hive can harvest in a day. Out of the blooming period, R bas < R max represents the basal resource available in the natural habitat during all spring and summer. When the fraction of natural habitat is small, a small increment of its value induces a proportional increase of the harvest. But when this fraction is large, the gain does not increase as much with the addition of natural habitat because it is limited by the foraging capacity of the bees. Let us say that x represents the fraction of natural habitat within the flight range of the hive and k s is the saturation constant of the harvest. The gain is:

G t = n f R max , during the crop blooming period, n f R bas , otherwise, (8)
where

R bas = R max x(1 + k s ) x + k s . (9)
The number of visits from the hive to the crop each day is proportional to the number of foragers and to the fraction of monoculture:

V t = n f (1 -x), (10)
where (1x) represents the fraction of crop within the flight range of the hive.

Results

Let us first explore the effect of competition on the pollination service. In previous models, the absence of competition produces a linear dependence of visits on the number of hives for a given map [START_REF] Hussler | Pollinator population size and pollination ecosystem service responses to enhancing floral and nesting resources[END_REF][START_REF] Lonsdorf | Modelling pollination services across agricultural landscapes[END_REF][START_REF] Olsson | Modeling pollinating bee visitation rates in heterogeneous landscapes from foraging theory[END_REF]. In our model the competition induces a saturation of the number of visits and even a decrease of the visits due to competition-induced hive mortality. We explore how the carrying capacity of . The carrying capacity of the resource is set as K = 7 g, and the map is the same as the one used in Fig. 1. The crop blooming period is from day 20 to day 55, with day 0 the first day of spring. The model ran 2 years with 100 hives placed at random on the border of the natural habitat, and we took the values for the second year. The corresponding plots for honey and for number of bees can be found in the Supplementary Material.

the resource and its renewal rate affects this dependence and thus the pollination efficiency regarding the number of hives placed in a crop on two maps (see Fig. 1). We can see that when the resource is in excess and the renewal rate is high, the number of visits is indeed proportional to the number of hives. When the resource is low, we can observe a saturation of the number of visits due to smaller hives, or even because some of the hives did not get through winter.

With this analysis we are able to estimate for each map an optimal number of hives for a given carrying capacity and renewal rate of the resources. We can see that the optimal number of hives increases with the amount of natural habitat, the carrying capacity and the renewal rate of the resource, but these effects are nonlinear and thus harder to predict.

In Fig. 2 we can see that the renewal rate of the resource does affect the number of visits only if it is really low, which can happen at the end of the blooming period for instance. We can conclude that most of the competition occurs within a day. The corresponding plots for honey production and number of bees can be found in the Supplementary Material, and they show a similar response.

We also explored the effect of the exponent γ of the site election function (Eq. ( 4)) on the mean foraging distance of bees and the mean gain of nectar and pollen per site. The results are shown in Fig. 3. We can see that the exponent performs as expected, increasing the foraging efficiency of the hives and reducing the foraging distance.

To quantify the effect of the amount of natural habitat we performed extensive simulations, monitoring the foraging activity of the bees on one thousand random maps during 4 years. For each map, we placed 100 nests at random on the borders of the natural habitat (which is the usual practice for managed hives). Due to runtime constraints, each map was simulated once. The number of maps was chosen so that the effect of the stochasticity in the heterogeneity of the landscape does not affect significantly the average results. The mean anual number of visits to the crop is shown in Fig. 4 for three possible blooming periods during the season. The corresponding contour plots for honey production and number of bees can be found in the Supplementary Material.

The number of visits to the crop reaches a maximum (darkest green in the plots) that remains stable for all the blooming periods (only its value increases when we go later in the season). We can also see that there is an additional dependence on the edge density of the natural habitat with a maximum between 0.2 and 0.4. A low edge density means that the natural habitat and the crop are well segregated, making some sites at the center of the crop unreachable for the bees. On the other hand, a large edge density means that the natural habitat and the crop are well mixed, diluting the resource of the natural habitat out of the blooming period. For the number of visits to be maximal, the patches of natural habitat can not be too dispersed, but cannot be extremely segregated either.

Results of the mean field model

As mentioned before, a mean field model is suitable to estimate the influence of the of natural habitat on the pollination service by A. mellifera without considering the effect of the distribution and shape adopted by this fraction. Here we present the corresponding results. We made simulations for different values of the fraction of natural habitat, and measured the number of visits of a hive during a 4 years period. The results are shown in Fig.

We can see that the production of honey increases with the fraction of natural habitat, saturating for large values of it. We can also see that the honey production barely depends on the blooming period of the monoculture. The number of visits, instead, strongly depends on the blooming period. It has a maximum for a fraction of natural habitat which ranges from 11.6% for an early blooming to 28.1% for a late one. The numbers of visits on late blooming crops are larger because the honey bee population is also larger during this time of the year.

Discussion

The spatially explicit model presented in this work provides a valuable insight about the impact of the landscape structure on the pollination service. We have seen that a mean field approach accurately predicts the impact of the fraction of natural habitat on the number of visits. Nevertheless, it fails to explain much of the richness of the possible outcomes since the disposition of the natural habitat patches and the competition between hives showed to have an influence on the pollination service. Even if the existence of an optimal fraction of natural habitat to maximize the number of visits is quite predictable with simpler mean field models or good sense, the value it takes in function of the configuration of the landscape, or the growth of hives remains quite non trivial. It is remarkable that some published experimental results suggest a maximum of pollination around the same values of edge density and fraction of natural habitat as the ones found by means of our model with default parameters (20% of semi-natural habitat for an edge density between 0.3 and 0.4) [START_REF] Martin | The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe[END_REF]. Moreover, this is, as far as we know, the only model which can effectively address the effect of the size of patches on the number of visits considering a cell size smaller than the size of foraging sites, allowing to find the explicit dependence on patch density. We could infer from extensive computations how this optimum varies with forager mortality or flight range of the bees, for instance. The model also showed that the moment of the flowering period does not influence much the optimal composition or configuration of the agricultural landscape for pollination.

Our model does not allow to predict an absolutely trustworthy number of visits. Indeed, several parameters (such as the number of trips per day and the number of flowers visited per travel), which have been considered constant, are in fact highly variable depending on parameters we chose not to investigate in this first study (weather, subspecies of bees or type of crop, for instance). But taking into account that these parameters vary evenly in every map, we can consider that the model is suitable to compare the number of visits between maps. Specific implementations for real world applications are under way, and will be the subject of further work.

The partial knowledge and the visibility of bees is hard to model in simple mathematical terms. That is why the site election modality of this foraging model, even if it appears less biased than in Bumble-BEEHAVE, remains an approximation. It clearly does not replace approaches based on individual decision making like the BEESCOUT to understand the foraging behavior of bees. But for the moment the computation time required to apply the decision making and field exploration to as much individuals as we have in a 3 km by 3 km map remains too high, and our function remains a good approximation. As it is shown in Fig. 3, the exponent γ can easily be deduced from the mean daily gain of foragers and the mean flight distance of the bees, and allows the user to avoid a quite unrealistic systematic choice of the most resourceful patch, leading to an overestimated competition between foragers for a single patch.

Our results apply to hives which are managed at low cost without intervention during spring and summer. In some places, with a particular lack of natural habitat offering flower diversity, agreements are made between farmers and beekeepers in order to keep alive enough hives for pollination purposes [START_REF] Goodwin | Pollination of crops in Australia and New Zealand[END_REF]. They can be moved across the fields through the seasons, and regularly supplied with syrup if needed. A good beekeeper can effectively erase the effect of landscape structure, but at a great cost in time and money.

In our analysis we have presented the total number of visits over the map as one of the important outputs, which is different from the economically more relevant fruit production. It is known that generally around ten visits per flower are enough for the fecundation to occur, and that more visits not only do not improve production, but can even damage the pistil and reduce fecundation probability [START_REF] Sáez | Extremely frequent bee visits increase pollen deposition but reduce drupelet set in raspberry[END_REF][START_REF] Rollin | Impacts of honey bee density on crop yield: A meta-analysis[END_REF]. We have also calculated production from the local number of visits each day, but because hives have identical foraging behavior, fruit production and number of visits are similar. Consequently, we did not find it useful to show these results here. However, we can say that the slight differences observed between visits and production only occur when the competition is strong and the flowers are overvisited. These differences are likely to increase for a larger number of hives. A detailed analysis will also be presented elsewhere.

It is important to rise a word of caution about our results concerning the impact of the fraction of natural habitat on production. What we demonstrated is that there is an optimal fraction of natural habitat needed in order to maintain a sufficient population of A. mellifera to pollinate the field. But it does not mean that increasing the fraction of natural habitat in an agricultural landscape will necessarily diminish the production. Indeed, a natural habitat such as a forest not only serves as a nectar supply for A. mellifera all over spring and summer, but also provides several other ecological services, e.g. as a wind breaker, air cooler, furnishing water retention and shelter for wild pollinators [START_REF] Marco | Services performed by the ecosystem: forest remnants influence agricultural cultures' pollination and production[END_REF], etc. The minimal fraction of natural habitat which begins to damage the pollination service by A. mellifera is likely to be lower than the one that will damage the actual yield of crop.

It is generally acknowledged that, while A. mellifera is the most used pollinator of crops, they are not the most efficient ones [START_REF] Garibaldi | Wild pollinators enhance fruit set of crops regardless of honey-bee abundance[END_REF][START_REF] Rollin | Impacts of honey bee density on crop yield: A meta-analysis[END_REF]. A model of pollination by wild bees using the same resource maps has already been developed by us, and will be used to characterize the effects of the landscape structure on the pollination by native pollinators on artificial and real landscapes, as well as the impact of mixing native and managed bees on pollination.

The goal reached by this work is a characterization of the effect of some structural features of landscapes on the pollination service. The model has been developed to analyze the pollination service in real GIS maps with a larger number of land covers and it is also able to predict the impact of slight changes in real landscape composition and configuration on the pollination service. It will be used in a subsequent study to predict the pollination service on different real agricultural landscapes. We hope that it can become a useful tool for farmers, beekeepers and policymakers in order to understand the impact of the composition and configuration of the agricultural landscape on the pollination service by the honey bee, help them to design their farms and contribute to a sustainable use of the environment.

Table 2: Summary of the main parameters used in the model. Other parameters have been set such that the results fit into a biologically relevant range of values for the number of bees and the production of honey: between 0 and 70,000 bees per hive, and between 0 and 50 kg of honey per hive, at the end of summer. Population dynamic parameters differ for each particular climate. We used a set of parameters that corresponds to the region of Río Negro (Argentina). Hive dynamics has been modeled following [START_REF] Khoury | Modelling food and population dynamics in honey bee colonies[END_REF], with parameters and a sensitivity analysis provided in the Supplementary Material.

Description

Figure 1 :

 1 Figure1: The effect of local competition. Left: A synthetic map, with 4% natural habitat (dark green) and edge density 0.19. Center: Visits to the crop. Right: Honey production. Each color corresponds to a value of the carrying capacity per cell, as shown. Resource renewal rate, r, is 100% per day. The crop blooming period is from day 20 to day 55, with day 0 the first day of spring. The model ran 2 years, and we took the values for the second year. The corresponding plot for number of bees can be found in the Supplementary Material.

Figure 2 :

 2 Figure2: Average number of visits as a function of the daily renewal rate of the resource (with minimum and maximum shown as error bars). The carrying capacity of the resource is set as K = 7 g, and the map is the same as the one used in Fig.1. The crop blooming period is from day 20 to day 55, with day 0 the first day of spring. The model ran 2 years with 100 hives placed at random on the border of the natural habitat, and we took the values for the second year. The corresponding plots for honey and for number of bees can be found in the Supplementary Material.

Figure 3 :

 3 Figure 3: Effect of the exponent γ of the site election function on foraging. The mean gain and mean distance are averaged over 10 realizations with a single hive placed at random, and the minimum and maximum values obtained are displayed as error bars around the mean. These simulations correspond to a map with only natural habitat, with carrying capacity of 7 g of resource and 100% resource renewal rate.

Figure 4 :

 4 Figure4: Mean annual number of visits to the crop as a function of the fraction of natural habitat and the edge density. The color scale shows the mean annual total number of visits. Left: crop blooming period from day 22 to day 55. Middle: crop blooming period from day 80 to day 115. Right: crop blooming period from day 130 to day 165. Day 0 is the first day of spring. We used 1000 maps for each of these analysis. Day 0 is the first day of spring. The corresponding plots for honey production and number of bees can be found in the Supplementary Material.

Figure 5 :

 5 Figure 5: Analysis of the impact of the fraction of natural habitat on the honey production and the annual number of visits with the mean field model. Top: Mean annual production of honey as a function of the fraction of natural habitat. Bottom: Annual number of visits to the crop as a function of the fraction of natural habitat. The colors correspond to different blooming periods, as shown in the legend in days. Day 0 corresponds to the first day of spring.

Table 1 :

 1 A summary of the main characteristics of several pollination-foraging models, as discussed in the Introduction.

		BEEHAVE	Olsson et al.	BEESCOUT	Hussler et al.	Bumble-	This work
			(2015)		(2017)	BEEHAVE
	Species studied	A. mellifera	Pollinators in	A. mellifera	B. terrestris	B. terrestris	A. mellifera
			general
	Level of	Colony	Community	Individual	Colony	Colony	Colony
	organization
	Spatial resolution	×	25 m / 3 km	Defined by user	25 m / 3 km	25 m / 5 km	10 m / 3 km
	Pollination feedback	×	×	×	×	×
	Inter-hive competition	×	×	×	×
	Site election	based on	Deterministic	Based on	Deterministic,	The best patch	Each patch has
	modality	detection	only based on	simulation of	based on	is chosen out of	a probability to
		probabilities	distance.	empirical flight	resource and	a	be chosen based
		derived from		patterns	distance.	probabilistically	on distance,
		BEESCOUT				chosen	resource and the
		and recruitment,				metapatch,	resource of
		which is affected				depending on	neighbor
		by the energetic				resource and	patches.
		efficiency of the				distance.
		trip
	Goal of the	Impact of many	Effect of	Understand	Effect of	Impact of many	Effect of
	model	stressors on	landscape on	landscape	landscape on	stressors on	landscape on
		individuals and	flower visit rates	exploration by	flower visit	individuals,	flower visit
		colony.	and bees' fitness.	bees.	rates.	colony,	rates.
						populations and
						community with
						mapping.

Acknowledgements

We acknowledge financial support several sources: FON-CYT, Plan Argentina Innovadora 2020: PICT 2016-0305, UNRN: PI 40-B-567, UNCuyo: SIIP 06/C546, CONICET: PIP 112-2017-0100008CO. We also thank the ENS de Lyon which allowed the visit of J. Joseph to Bariloche through their internship program. We also thank the reviewers of the manuscript for their valuable comments and suggestions.

The funding agencies were not involved in the research or the preparation of the manuscript.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: