Quoc Tung Lam Nguyen

Tuan Tran
email: quoctuan.tran@cea.fr

Yvon Besanger
email: yvon.besanger@g2elab.grenoble-inp.fr

Hoa Van

Nguyen

Van Hoa Nguyen
email: vanhoa.nguyen@cea.fr

Tung Lam Nguyen
email: tung-lam.nguyen@g2elab.grenoble-inp.fr

Synchronization conditions and Real-time constraints in co-simulation and Hardware-in-the-Loop techniques for Cyber-Physical Energy System assessment

Keywords: Co-simulation, Hardware-in-the-loop, Synchronization Condition, Real-Time constraint, Time semantic

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

In the course of modernization the power system towards a decarbonization scenario, it is required to have a high penetration of Distributed Renewable Energy Sources (DRES) and to develop more advanced Information and Communication Technologies (ICT) and automation solutions, including big data assessment [START_REF]Smart grids in distribution networkroadmap development and implementation[END_REF]. The integration of DRES gives rise to new challenges in maintaining the security of supply and the power quality due to their intermittent behaviours. It leads to new rigorous requirements in testing and validating innovative solutions in the integration of such devices into the power systems [START_REF] Jong | European white book on real-time power hardware in the loop testing[END_REF]. It is even more complicated because the power system itself has become a sophisticated Cyber-Physical Energy System (CPES) with complex inter-connection and inter-dependency. Sole component testing outside of its deployment environment (power, ICT, automation, etc.) may not allow us to fully assess the impact of DRES to the CPES. A holistic validation approach is necessary to ensure power system stability assessment and rapid prototyping and to validate safety and security of power system operation with innovative solutions.

Achieving such a holistic approach is however not an easy task due to several important challenges:

• Firstly, the existing domain-oriented simulation tools and prototypes fail to address the complex behaviour and the inter-dependency of the multi-domain CPES (Electricity, Thermodynamics, ICT, market, etc.).

• Secondly, the CPES contains applications with various time scales (e.g. transient and steady state operations) which is challenging to be included in one single experiment with existing simulation tool without overloading the simulator.

• The consideration of the interaction of these elements requires an appropriate method that includes coupling of continuous (energy behaviours) and discrete events (ICT system or market commands).

• Inaccurate or incomplete models that is preferred to be replaced by their real-counterpart in an experiment. This approach is however not very flexible when we have to consider faulty and extreme conditions.

To cope with these challenges, several innovative approaches for assessment and validation of CPES can be considered: cosimulation, Real-Time Simulation (RTS), Hardware-in-the-Loop (HIL) and the combination of those techniques. Due to the multi-domain nature of smart grid, these new techniques involve the combination of several entities (simulation, physical component, communication network, etc.). An important aspect in the implementation of these technologies is the coupling of different simulations or processes and the coupling of hardware-software. Besides the necessity of harmonization of different data formats and protocols, in order to correctly reflect the real scenario, it is mandatory to ensure that all the elements of the framework share the same time horizon and semantics [START_REF] Palensky | Cosimulation of intelligent power systems: Fundamentals, software architecture, numerics and coupling[END_REF][START_REF] Nguyen | Using power-hardware-in-the-loop experiments together with co-simulation in a holistic approach for cyber physical energy system validation[END_REF]. The constraint is even stricter when a physical equipment is involved in the test, because the "virtual" elements need to keep up with the real-time behaviour, which is particularly challenging due to the capacity of the simulator and due to the potential variability of time steps employed by the solvers [START_REF] Nguyen | Using power-hardware-in-the-loop experiments together with co-simulation in a holistic approach for cyber physical energy system validation[END_REF].

The synchronization among the individual entities, so that they can represent correctly the desired scenario (the real deploying environment or an extreme situation, for example) was not sufficiently considered in the literature. In actual literature, the coherency of the processes in a validation framework is ensured by applying coupling and synchronization techniques (cf. section 3) establishing arrangement of events among the participating entities. The existing techniques synchronize the processes based on time steps while the meaning and the semantic of time among the simulators/hardware are manually managed by the system designer.

It is therefore important for the engineers to have a theoretical formalization of such synchronizations which includes timing semantic [START_REF] Chang | Heterogeneous simulation -mixing discrete-event models with dataflow[END_REF][START_REF] Cremona | Hybrid co-simulation: It's about time[END_REF][START_REF] Lee | What is real time computing? a personal view[END_REF]. On one hand, it helps to determine whether a framework is feasible, allows configuring the experiment systematically to enhance and ensure timing consistency. On the other hand, it provides a theoretical base to extend the synchronization conditions to more complex situations, e.g. variable step simulation or integration of HIL to co-simulation in presence of latency. In this paper, we attempt to fill in this gap with an original formulation of consistency condition that can be used to build and to verify the semantics of such a validation framework. The condition detaches itself from the test implementation techniques and focuses on finding the closest to reality synchronization among the entities. From these results, we revisit the real-time (RT) constraints in the integration of HIL to co-simulation framework and propose a relaxed condition that can be used to configure a HIL framework and to seek for an optimal compromise between the time resolution and the synchronization strictness.

The paper is organized as following: in section 2, we provide a brief review on the advanced technologies in the validation of CPES, followed by a certain focus on their associated synchronization techniques in section 3. In section 4, we formulate the problem with time semantic in a coupling interface and the consistency conditions (synchronization and real-time constraints) are then proposed. The application of the conditions is then illustrated via a remote HIL test-case for secondary control in islanded microgrid with multi-agent system in section 5. The paper is then concluded in section 6.

Advanced CPES validation technologies

CPES, as a system of systems with diverse natures, requires the validation framework to reflect the integrity and the inter-dependency of the subsystems. A single domain solution may not be the correct solution for the whole framework. We review in this sections advanced CPES validation technologies : cosimulation, RTS and HIL as well as the approach of laboratory coupling.

Co-simulation.

A co-simulation is composed of coupled simulators, running on their own solvers and models, that cooperate with each other. The output of one simulator is coupled to the input of another and vice versa for dynamic exchanges. Cosimulation is used in CPES for coupling of distributed model in the same domain (e.g. real-time ElectroMagnetic Transients (EMT) simulation of smart grids [START_REF] Guo | Comprehensive real-time simulation of the smart grid[END_REF], hybrid (SPICE-type) circuit EMT stability simulations [START_REF] Asghari | Intefacing techniques for electromagnetic field and circuit simulation programs[END_REF], and paralleled EMT simulations [START_REF] Zhou | Parallel massive-thread electromagnetic transient simulation on gpu[END_REF]) or between sub-systems of different natures (e.g. analysis of wide area monitoring and control [START_REF] Lin | A study of communication and power system infrastructure interdependence on pmu-based wide area monitoring and protection[END_REF], control and optimization in distribution networks [START_REF] Bottura | Ict-power cosimulation platform for the analysis of communication-based volt/var optimization in distribution feeders[END_REF][START_REF] Armendariz | A co-simulation platform for medium/low voltage monitoring and control applications[END_REF][START_REF] Fennibay | A heterogeneous simulation and modeling framework for automation systems[END_REF][START_REF] Zhabelova | Multi-agent smart grid automation architecture based on iec 61850/61499 intelligent logical nodes[END_REF][START_REF] Wei | Greenbench: A benchmark for observing power grid vulnerability under data-centric threats[END_REF], and distributed energy integration [START_REF] Dufour | On the use of real-time simulation technology in smart grid research and development[END_REF][START_REF] Burnier De Castro | Dynamic co-simulation of agent-based controlled charging electric vehicles and their impacts on low-voltage networks[END_REF], impact of communication latency on the power system [START_REF] Stifter | Steady-state cosimulation with PowerFactory[END_REF][START_REF] Karmi | A flexible smart grid co-simulation environment for cyber-physical interdependence analysis[END_REF] and impact of cyberattacks on the electric power grid [START_REF] Wei | Greenbench: A benchmark for observing power grid vulnerability under data-centric threats[END_REF][START_REF] Mallouhi | A testbed for analyzing security of scada control systems (tasscs)[END_REF]).

Real-time Simulation. RTS in CPES is the concept of reproduction of the output voltage/current waveforms representing the behavior of the real power system being modelled, with desired accuracy in term of timing (with respect to physical time), amplitude and frequency, etc [START_REF] Lauss | Characteristics and design of power hardware-in-the-loop simulations for electrical power systems[END_REF]. In the context of CPES assessment, RTS is important for the validation of timing-sensitive safety-critical applications [START_REF]Real-time simulation technologies for power systems design, testing, and analysis[END_REF][START_REF] Dufour | On the use of real-time simulation technology in smart grid research and development[END_REF], usually in HIL setup.

Hardware-in-the-Loop. (HIL) technique can be considered as a hybrid configuration with software (simulation) and hardware (physically connected) interfacing through digital and analog input/output (I/O) signals. HIL in CPES is widely used for testing of single components (e.g. protection relays [START_REF] Almas | Over-current relay model implementation for real time simulation & Hardware-in-the-Loop (HIL) validation[END_REF], power system devices [START_REF] Loddick | The use of real time digital simulation and hardware in the loop to de-risk novel control algorithms[END_REF], novel power converter structures [START_REF] Kotsampopoulos | Design, development and operation of a PHIL environment for Distributed Energy Resources[END_REF]) or for CPES as a whole [START_REF] Bian | Realtime co-simulation platform using OPAL-RT and OPNET for analyzing smart grid performance[END_REF][START_REF] Nguyen | Agent based distributed control of islanded microgrid real-time cyberphysical implementation[END_REF].

RTS and HIL can also be combined together to create a complete and realistic environment for CPES validation [START_REF] Nguyen | Using power-hardware-in-the-loop experiments together with co-simulation in a holistic approach for cyber physical energy system validation[END_REF][START_REF] Nguyen | Real-time simulation and hardware-in-theloop approaches for integrating renewable energy sources into smart grids: Challenges & actions[END_REF]. Laboratory coupling. CPES is complex system of systems that requires adequately sophisticated validation framework that sometimes exceeds the capacity of a single RI, in terms of equipment and expertise. Combining complementary expertise and platforms of different Research Infrastructures (RI) in a holistic test is therefore a judicious solution that does not require heavy investments in time, infrastructure and human resources. Implementing cross-infrastructure holistic testing allows the RI to exchange specialized expertise and to get access to a shared resource pool, to manually or remotely borrow or share the infrastructure for research activities. Besides, this approach provides also the technical base and paving the way to international collaboration by combining several infrastructures and/or replacing non-available components/systems by simulation, increas-ing the realism of innovative validation methods [START_REF] Nguyen | Using power-hardware-in-the-loop experiments together with co-simulation in a holistic approach for cyber physical energy system validation[END_REF][START_REF] Lehfuss | Coupling of real-time and co-simulation for the evaluation of the large scale integration of electric vehicles into intelligent power systems[END_REF][START_REF] Stevic | Multi-site european framework for real-time co-simulation of power systems[END_REF][START_REF] Lundstrom | Evaluation of system-integrated smart grid devices using software-and hardware-in-the-loop[END_REF][START_REF] Marten | Result processing approaches for large smart grid co-simulations[END_REF]; from local scale to multi-platform scale.

Synchronization techniques and time semantic

The implementations of co-simulation, parallel computing, hybrid simulation and hardware-in-the-loop depends on synchronization techniques to ensure the correctness and consistency of the validation framework with respect to the real scenario. Existing works in the literature focus on the technical coupling aspects, in term of formulating master algorithm and building interfaces. Two main approaches can be acknowledged: conservative approaches (e.g. Null Message algorithm [START_REF] Chandy | Distributed simulation: A case study in design and verification of distributed program[END_REF] or Global Synchronization approaches, such as: Bounded Lag algorithm [START_REF] Lubachevsky | An analysis of rollbackbased simulation[END_REF], Time Buckets [START_REF] Steinman | Breathing time warp[END_REF] or Composite Synchronization [START_REF] Nicol | Composite synchronization in parallel discreteevent simulation[END_REF]) where the concurrent processes can only advance for a time period in which it can be proven that no events sent out from others are expected (i.e. "lock step" approach); and optimistic approaches where sub-processes are allowed to advance concurrently, but backward regression must be performed once a conflict is registered [START_REF] Jefferson | Virtual time[END_REF]. Actual implementations of conservative approaches depend on the nature of the models (discrete event or continuous) and the way sub-processes advance their time steps (uniformly or variably). For coupling of continuous simulations, Gauss-Seidel (or Serial Sequence) [START_REF] White | Waveform relaxation: Theory and practice[END_REF] and Jacobi (or Parallel Sequence) [START_REF] Sand | A jacobi waveform relaxation method for ode[END_REF] techniques can be applied [START_REF] Palensky | Cosimulation of intelligent power systems: Fundamentals, software architecture, numerics and coupling[END_REF]. Master-Slave, time-stepped, global event-event driven are popular synchronization technique for coupling of discrete event (DE) simulations or DE -continuous couplings [START_REF] Weilin | Co-simulation for smart grid communications[END_REF][START_REF] Nguyen | On conceptual structuration and coupling methods of co-simulation frameworks in cyber-physical energy system validation[END_REF]. Another interesting approach for strong coupling of simulators, especially in the presence of communication latency, is waveform relaxation method (WRM) [START_REF] Nguyen | On conceptual structuration and coupling methods of co-simulation frameworks in cyber-physical energy system validation[END_REF], where an individual simulator is represented by complete I/O waveform at the coupling points, to the other simulators. Based on this assumption, the simulation is solved independently inside the two simulators and the results are exchanged to be input for a following iteration. The process stops when the convergence conditions are reached. WRM is generally slower than others due to convergence time. However, when there is a non-neglectable latency among the individual entities, WRM can reach the final result much faster than the other methods because it requires much fewer number of exchanges.

The existing synchronization methods can be considered as arrangement techniques that help to chain the events among simulators to a desired order. The meaning of that order (i.e. the desired time horizon of the system) is (implicitly or explicitly) defined and ensured by the user. Besides, they actualize synchronization based on the basic computational cycle of a simulator. The synchronization step can be either micro-step or marco step, but the step is often supposed to be fixed or predictable to facilitate the synchronization. Moreover, many of the existing techniques require mechanisms to stop or roll back individual simulator to execute the synchronization. This is however not always allowed in case of asynchronous cosimulation or interaction with hardware/real environment.

Formulation of consistency conditions

Problem Formulation

To facilitate the reading, the following conventions are used in the following:

• A normal character represents a specific variable (e.g. t i represents a particular moment).

• A capital character represents the set of variables (e.g. T represents the set of t i ).

• A bold character represents a mapping or a function (e.g. t M ).

Time step (t s ) is defined by the model or the user to capture all the fluctuations of the simulated system and therefore can be constant, fixed or variable depending on the choice of model of computation and solver. Time step in CPES simulation can range from microseconds to minutes [START_REF] Mets | Combining power and communication network simulation for cost-effective smart grid analysis[END_REF]. To accurately consider transient behaviours of a 50-/60-Hz power system, a time step of roughly 50µs or smaller is required. On the other hand, the system level interaction may allow longer time steps [START_REF]Real-time simulation technologies for power systems design, testing, and analysis[END_REF]. The continuity and variability of time steps does not indicate if a model is discrete or continuous, but rather the manner which their models of computation handle time. e.g. a discrete time continuous model can have a variable time step for specific solver usage, or fixed step according to user's request (e.g. in Matlab), whereas, in discrete event simulation, time advances only when an event occurs, uniformly or discretely. Definition 1. The simulation horizon of a model is defined as the subset of physical time (i.e. wall-clock or world time) that is simulated/represented by the said model.

T M = {t ∈ T W t ≥ inf t m ∧ t ≤ sup t m }
In other word, it is the part of physical time "represented" by the set of the model time steps.

The relation between T M , t s , τ is decided by the simulator according to a particular need (i.e. accelerate the simulation, get more details from a specific phase, adapt to the solver or actualize real-time simulation and hardware-in-the-loop, etc.).

We characterize the mapping from τ to T M of a model by a function t M .

t M ∶ τ M → T M τ n ↦ t M (τ n ) = ∑ {i∈τM i<n} T s (τ i ) (1)
The function t M characterizes the way a simulated model advances its simulation horizon, not necessarily in uniform or linear with world time (cf. figure 1). t M is linear if the time step t s is constant and fixed and is non-linear in case of variable t s . Definition 2. C M is defined as the set of time stamps in T W on which the simulator advance from one computational cycle to its following computational cycle. The mapping from τ M to C M is characterized by the function:

C M ∶ τ M → C M τ n ↦ C M (τ n ) = ∑ {i∈τM i<n} T E (τ i ) (2)
C M (τ n ) represents the elapsed (world) time that the simulator has spent to reach the tag τ n . In general, T E (τ i ) (and as a consequence C Mi ) can not be exactly defined as it varies depending on the complexity of the computational cycle, the efficiency of the solver/simulator and the power of the computer (cf. figure 1). However, T E can be observed and estimated to some degree of accuracy, for a certain set of model, solver and simulator, with respect to a specific computer.

In general, real-time capacity of a CPES simulator is indicated following the condition 1.

Condition 1. (Necessary Condition) A simulator is considered to be Real-Time compliant with a particular model if the Execution time of that model on the said simulator is less or equal to the its time step:

∀τ i ∈ τ M , T E (τ i ) ≤ T s (τ i ) (Sufficient Condition) A simulation is Real-Time if its computa- tional cycle,
including idle time is equal to its time step:

∀τ i ∈ τ M , T E (τ i ) + T idle (τ i ) = T s (τ i )
where the set of tags τ is defined as the subset of N that indexes the sequences of time steps, with which the model is simulated, T E is the World Time required to finish the τ i -th computation cycle; T s is the incremental interval in which the state of the system is simulated and advanced. In a simpler point of view, τ can be considered as the set determining the order of t s .

The set τ defines a semantics that allows the simulator to reason about the ordering relationships and to detect inconsistencies. Two events are "synchronous" if they have the same tag (i.e. happen at the same time step), or else, one unambiguously precedes the other. We can denote the t s associated to the tag τ i as t s (τ i ), which signifies that t s can also be variable. A simulation is considered Offline when the execution time differs from the selected time step. If T E < T s , the simulation is classified as Fast. Real-time capacity is enabled by putting the simulator Idle until the next time step. On the other hand, if T E > T s , overruns occur and the simulation is classified as Slow. In that case, either the model has to be simplified (not always possible) or the time step has to be increased (less information/details) to enable the Real-time capacity.

Consistency conditions 4.2.1. Synchronization conditions

The synchronization techniques ensures the correct mapping of two simulators. However, they synchronize the computational cycles and not the simulation horizons (i.e. based on the assumption that their simulation horizons and the ways that they advance states are coherent.). In order to have a correct temporal semantic for a co-simulation framework, one needs not only the correct mapping of the steps, but also the coherency between the time steps and the simulation horizons, especially in case of variable time steps or when then simulation horizon does not advance uniformly. We envisage in this section, a generic condition of synchronization that take into account all these elements. The proposed condition, will be a useful tool for the user to design and synchronize a co-simulation framework.

Considering a co-simulation interface between simulators

A and B, we propose also the following definitions: Definition 3. S AB is defined as the set of points in T W that A updates its states with respect to B. R AB is defined as the set of points in T W that B reads the states of A.

S AB and R AB represent the rate with which A exhibits its internal state to B and the rate with which B update its internal variables (copies) corresponding to the state of A (cf. figure 1). In general, S AB and R AB are often associated to the I/O task of the micro-step or marco-step and are often chosen as same points for all the co-simulation interfaces (i.e. in case of ad-hoc co-simulation of several simulators). It is however necessary to distinct them in the formulation because they are actually different tasks and can be separated accordingly (i.e. instead of updating the state every cycle, A can maintain the state internal and only update once every several steps; same for B, the simulator can choose to read the state of A once in a while instead of doing it every cycle). Besides, these definitions englobe both cases of synchronization with micro-step or marco-step.

We characterize also the mapping functions from the computational timelines (C A , C B ) to S AB and R AB .

S AB ∶ S AB → C A s ab ↦ S AB (s ab ) = max {ca∈CA ca<sab} c a (3)
which means that the updated state of A at s ABi point is the latest calculated state of A before s ABi .

R AB ∶ R AB → C B r ab ↦ R AB (r ab ) = min {cb∈CB cb>rab} c b (4)
which means that the state of A received at point r ABi will be acknowledged in the very next cycle of B following r ABi .

In a similar manner, we consider the mapping from S AB to R AB . Certainly, it is possible that B chooses a longer (faster) refresh than the rate A exhibits its state to B. In this case, the copy of state of A in B will be updated slower (faster) and there is information loss in term of state of A (or the state will be duplicated). This has to include also the communication time between the two simulators, if that one exists.

M AB ∶ R AB → S AB r ab ↦ M AB (r ab ) = max {sab+TcomAB<rab} s ab (5)
An illustration of these functions in a co-simulation interface from A to B is given on the figure 1. In general, S AB ,R AB and M AB characterize the exchange from A to B. They help to determine the actual computational cycle that we get/set the state of A. In a similar manner, we can define the functions characterizing the exchange from B to A: S BA ,R BA and M BA . It is important to note that these are not the inverse functions of S AB ,R AB and M AB as the signal flows can be configured differently from A to B and vice-versa.

Their inverse functions, on the other hand, allow us to determine the points where the forward functions get values from.

S -1 AB ∶ C A → S AB c a ↦ S -1 AB (c a ) = min {sab∈S AB sab>ca} s ab (6) R -1 AB ∶ R AB → C B c b ↦ R -1 AB (c b ) = max {rab∈RAB rab<cb} r ab (7)
In the sense that if r ab is not updated yet, next cycle of B will still use the old value it received from A at precedent read point.

M -1 AB ∶ S AB → R AB s ab ↦ M -1 AB (s ab ) = min {rab-TcomAB>sab} r ab (8)
From these definitions, two things can be acknowledged: firstly, the formulation of S, R and M are for a general case and they may be specified in a slightly different way according to the operator's specifications (e.g. ignoring the latest state and gets the previous state). Secondly, the mapping in forward and inverse directions are not always bijective. Therefore, upon synchronization, it is possible that information loss may take place (i.e. we can fully reconstruct the time horizon of A from the information received at B only under some particular conditions. The formulation of these "loss-less" condition in cosimulation is also a very interesting subject but is not in the scope of this paper.

From the above formulation, we try to establish a generic condition for synchronization (i.e. on which condition we can verify that a synchronization is correct.). In order to know if a co-simulation is correct, we need to get back to the purpose of a co-simulation framework. The uttermost objective of such a framework is to combine the behaviors of several sub-systems to deduce the behaviour of the holistic system in a certain time horizon. Certainly, for the combination to make sense, it is required that the simulation horizons of these sub-systems are coherent.

I S ≠ ∅ (11)
For example, it is not right to integrate a simulation of a building consumption during the night with a simulation of PV energy production during the day, because their simulation horizon does not intersect and combining them will introduce an unrealistic scenario. However, if we include a battery simulation that is charged during the day by the PV and discharge during the night to supply the household, then the co-simulation framework is well-posed, because the simulation horizon of the battery will cover the disjoint part of PV and household (i.e. I S ≠ ∅). In the condition 2, it is possible that some simulators in S may start sooner or extend longer than I S , however in that case, these extensions may not include some other simulators and their interactions, i.e. the results may not be complete and accurate.

The condition 2 is a necessary condition, it does not ensure that the result of the co-simulation framework is accurate, i.e. the PV production at 9 a.m maybe mapped to charge the battery at 9.15 a.m. In order to ensure that the simulators get and send correct values at exact moments and use them in a proper manner to produce accurate result, we need a sufficient condition of synchronization. We propose in the following the generic condition under which the co-simulation is accurate.

As aforementioned, in digital simulation, engineers adopt the approach with discrete time, i.e. the state of a model is considered as unchanged during the time step. Therefore, considering a co-simulation interface between two models A and B with zero-order hold extrapolation, if we want to map a signal from A at t A to B at t Bi , it is necessary that the mapping from A to B of t A will land inbetween t Bi and t Bi+1 .

From the equations (1),(2)(3)(4)(5), we can define the mapping function from A to B (as illustrated on figure 1):

f A→B (t A ) = t B ○ C -1 B ○ R AB ○ M -1 AB ○ S -1 AB ○ C A ○ t -1 A (t A ) (12) It is worthy to take note that f A→B ≠ f -1
B→A and f A→B is not necessarily bijective.

In the ideal case, we want to map t A to the exact t A moment in T B , however, as T B is discrete, we need to map t A to the time step that contains t A in T B , i.e:

t Bi ≤ t Ai < t Bi+1 ⇒ t Bi ≤ f A→B (t Ai ) < t Bi+1 ( 13
)
Condition 3. (One-sided Condition) Model A is synchronized to model B, noted as A ↷ B iff ∀t A ∈ T A , if t Bi ≤ t A < t Bi+1 then t Bi ≤ f A→B (t A ) < t Bi+1 (14)
Similarly, as f A→B is not always bijective, if A ↷ B, then it is not necessarily true that B ↷ A.

A B ⇐⇒ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ A ↷ B B ↷ A (15)
The condition 4 requires that the interface between A and B satisfies the condition 3 in two directions, i.e.

∀t A ∈ T A , t Bi ≤ t A < t Bi+1 → t Bi ≤ f A→B (t A ) < t Bi+1 and ∀t B ∈ T B , t Ai ≤ t B < t Ai+1 → t Ai ≤ f B→A (t B ) < t Ai+1 .
It means that the mappings lead to correct time steps and does not mean necessarily T A = T B (i.e. not a bijective relation). This relaxation is particularly useful when we deal with simulators with variable time steps and is more suitable to be used in co-simulation framework with Jacobi method.

The condition 3 can be enforced by requiring also that the responses from the opposite model to return safely in one step. It leads to a stricter version of one-sided synchronization condition.

Condition 5. (Strict Condition) Model

A is synchronized to model B, iff ∀t Ai ∈ T A t Bi ≤ t Ai < t Bi+1 → ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ t Bi ≤ f A→B (t A ) < t Bi+1 t Ai < f B→A (t Bi+1 ) < t Ai+1 (16)
It leads also to a stricter version of two-sided synchronization. If we want to impose this latter (and stronger) restriction, we would have the classical "lock step" synchronization condition in literature that can be ensured by a Gauss-Seidel method. This condition is suitable to be applied to Co-simulation with Master-Slave scheme.

From the above analysis, we can come up with a sufficient condition as to the accuracy of the co-simulation framework as: Condition 6. (Sufficient Condition) Given a co-simulation framework S of n models S 1 , S 2 , . . . , S n , S is correctly synchronized iff ∀i, j = 1..n i ≠ j, S i S j

The conditions 2 and 6 constitute the requirements for a cosimulation framework to be correctly configured in a generic case. This allows us in one hand, to justify and to adapt the classical (micro-marco)step-based methods such as Gauss-Seidel and Jacobi in a particular test-case, and on the other hand, to extend them towards the possibility to work with variable timesteps or multi-time-scale applications. Moreover, as the formalization treats the solvers in their generality, it is possible to use these conditions with discrete event simulators where it is particularly difficult with classical approach due to its unpredictable steps.

For instance, as of for now, a roll-back mechanism is required when a continuous model is synchronized with a discrete event model, i.e. once the signal from the DE model is detected, if it does not correspond to the current step of the continuous model, then the model is rolled back and the predictable step is then adjusted to match the signal timing. It is particularly problematic to actualize, is often not supported by various softwares and will slow down the co-simulation significantly. As we relaxed the matching from one point to an interval, it would reduce the number of roll-back and improve the simulation performance. It can be eventually be possible to use the formulation to estimate and to tune the communication points to minimize the chance of roll-back.

We acknowledge also that allowing the matching from one point to an interval will also benefit the simulator's performance. As computer use floating number to model "time". It is particularly more difficult to detect the matching of time stamps than to compare them. Certainly, it introduces another uncertainty to a context that already dominated by uncertainties, but it provides us with more possibilities, such as synchronizing with variable marco-step. Moreover, a proper formalization allows us to relax the classical Real-time constraint to a more flexible condition that can be used for slow-simulation. We consider in next section, the formulation of a relaxed real-time constraint.

Real-time constraints

HIL experiment can be imagined as a "co-simulation" between a real-time simulator A and the real hardware. The difference is that the real-hardware "simulator" B is assumed to have infinite bandwidth, infinitesimal resolution and instantaneous computation/reaction time, i.e. T B ≡ T W , C B ≡ T W and t B , C B are bijective. S AB and R AB is however limited to the output frequency and the acquisition rate of the A/D converters or the sensors of the hardware. M AB can be either instantaneous or the latency (communication and power reaction) between the simulator and the hardware, depending on the decision whether the linkage is a part of the hardware under test (HUT) (e.g. does it reflect the real setup).

Basically the classical RT condition 1 requires that the loop time is safely bounded in one time step of the RT simulator. It can be considered as a special case of the strict one-sided synchronization condition 5 where B has infinitesimal resolution, i.e. B responses instantaneously instead of next time step. This condition is ensured inside the DRTS but as soon as the signal is emitted from the DRTS and enters the loop (communication to the hardware, response time and feedback to the DRTS), it is very hard to maintain the determinism, i.e. a feedback from one time step can be applied to a another time step later on, which is not correct. Until now, this in-exactitude is accepted as a part of the HIL uncertainty.

We can apply the conditions 2 and 6 to determine whether a loop is "real-time compliant" and see to what extent we can specify the parameters of the loop to satisfy it. The idea is relaxing the strict RT condition to a more flexible one. Instead of requiring the "time loop" to be bounded in a time step, we can use the condition 3 to ensure that received value is from the exact sampling period, not from previous period.

S HuT ⇐⇒ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ T S ∩ T HuT ≠ ∅ S ↷ HuT HuT ↷ S (18)
In the case that S is a DTRS, then T S ≡ C S . The condition can be simplified:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ T S ∩ T HuT ≠ ∅ ∀s S -HuT i ∈ S S -HuT , s S -HuT i ≤ M -1 S-HuT (s S -HuT i ) < s S -HuT i+1 ∀r HuT -S i ∈ R HuT -S , r HuT -S i ≤ M -1 HuT-S (r HuT -S i ) < r HuT -S i+1 (19)
It is particularly interesting. While we do not have much control over M S-HuT or M HuT-ST , S S -HuT and R HuT -S are the factors that we can configure and manipulate. The principle idea of the conditions 8 and 7 is the transformation from the strict conditions (5 and 6) to more simplified and practical conditions that can be easily applied in real set up (i.e. consisting of only variable that the users can freely manipulate). The simplification is quite straight forward. In the case of RT simulation, the functions t and C are ensured by the internal interrupters and the FPGA cards, so we only need to handle S, R and M. S and R can be manipulated by the user, while M can be estimated and bounded to a certain degree of exactitude (as in the following case-study). From the analysis on M S-HuT and M HuT-S , we can see to what extent we can correctly specify S S -HuT and R HuT -S .

Example: Cross-infrastructure secondary control of microgrid

We consider in this section a test case of Distributed Multi-Agent System (MAS) control of islanded microgrid. The testcase is implemented by integrating control agents in PRISMES platform (Le Bourget du Lac, France) with the RT Simulator (OPAL-RT) in PREDIS platform, Grenoble INP (Grenoble, France) -70 km apart.

System configuration

The distributed control scheme of microgrid based on MAS is previously proposed in [START_REF] Nguyen | Multi-agent system with plug and play feature for distributed sec-ondary control in microgridcontroller and power hardware-in-the-loop implementation[END_REF] and is composed of three layers (figure 2). The main differences of this novel architecture with the classical centralized approach are:

• An agent layer is added on top of the local control layer.

The role of agent layer is to actualize optimization and to bring the microgrid state to desired conditions (e.g. regulation of voltage or frequency). The agents affect the grid via sending set-point to the local controllers.

• An agent only communicates with its neighbourhood and does not necessarily require knowing the state of the whole microgrid. The stability and consistency of the grid is ensured via a consensus algorithm among the agents. Due to this consensus process, the agents only send set-points to the primary controller once their outputs converge.

• The agent is developed in a plug and play manner to promote interoperability. Once an agent is disconnected or a new one joins the network, the algorithms in other agents in its neighbourhood will detect automatically and adapt themselves in the next iteration. The considered experiment involves the integration of remote CHIL. The test-case is implemented with the virtualized MAS in PRISMES (CEA, France) and the OPAL-RT simulator in PREDIS, Grenoble INP (70). The considered network is given on Figure 3. The microgrid is simulated in real-time with the RT Simulator OPAL-RT in PREDIS platform. A Redis server is also hosted in PRISMES serving as real time exchange point. In this experiment setup, the agents are virtualized into python processes and are launched at CEA (France) (figure 4).

The following scenario was considered: Load 1 connects to the microgrid at 0s, load 2 connects to the microgrid at 60s and load 1 is disconnected from the microgrid at 120s.

Consistency condition analysis

In order to correctly configure the implementation, we consider the simplified real-time condition 7 with the real-time simulator OPAL RT (PREDIS) and the MAS (PRISMES) as the HuT. The synchronization is actualized via Master-Slave manner. While the communication between the MAS and Redis server is actualized via the local area network and benefits a consistent 4ms latency, the inter-platform communication requires a latency assessment (figure 5 and6).This latency characterizes the communication channel between the Redis server and the RT simulator (i.e. the M OPAL-Redis mapping). The histogram of the measured latency (minimum latency: 28 ms and maximum: 85 ms) is then constructed and is fitted to the basic Kernel distribution (figure 5) and the more biased Gamma distribution (figure 6). In the case with Kernel distribution fitting, the cumulative distribution function (CDF) shows us that in 99% of the cases, the latency is smaller than 62 ms (figure 5). The same analysis with CDF in the case with Gamma fitting gives the boundary of 99% confidence interval at 51 ms (figure 6). The difference demonstrates that the choice of distribution model plays an important role on the results of the analysis. Until now, the most appropriate distribution model to describe communication latency is still a controversial topic. Most of the time, people employ a simple interval or a normal distribution. The choice of the distribution models however affects the estimated bounds of the M function in the synchronization condition.

Considering that the data on Redis server is updated instantly (in fact there is a delay but it is neglectable with respect to the latency), S OPAL-Redis is characterized by the output rate of the RT simulator to the server and R Redis-OPAL is characterized by the request/reception rate of the simulator (round trip). The condition 7 is rewritten as:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ s OPAL-Redis (i) ≤ M -1 OPAL-Redis (i) < s OPAL-Redis (i + 1) r Redis-OPAL (i) ≤ M -1 Redis-OPAL (i) < r Redis-OPAL (i + 1) (20)
Since we analyzed that M -1 Redis-OPAL (i) and M -1 OPAL-Redis (i) can be bounded with 99% confidence at 62 ms, the advancement of S OPAL-Redis and R Redis-OPAL can be configured to make sure that the condition 7 is satisfied. Particularly, output step of the real-time simulator should be larger than 62 ms and the request rate should be larger than 124 ms (e.g. a rate of 150ms will work). Taking into account the other delay (MAS-Redis, Redis processing time and MAS processing time) and considering that secondary control does not strictly require very fast reaction, in the real setup, we chose to set both values to 500 ms. In this setup, the MAS system is "soft" synchronized to the RT simulator, which is set to run at much shorter time step (50 µs) to capture transient behaviour of the grid.

In the eventual case that the simulation loses synchronization due to high latency, it can be considered similar to an overrun. The following is totally a choice of the user, either stopping the experiment or continueing and compensating at the next sync point, if the small in-exactitude is tolerated. Particularly, in a private network and without multiple intermediate servers, the latency should be very stable.

Comparing to the hard synchronization approach, which requires much more strict setup, such different time steps with "soft" synchronization condition helps to avoid unnecessary computational load on one side and information loss on the other side. It may however cause potential desync while the framework is executed for a long period if ones time step is not a multiplier of the other; as the error may accumulate. Moreover, if the simulation are not real-time, a more difficult problem for synchronization is the difference in processing power, especially in case that an individual entity in the framework does not support stop or roll back.

Experimental Results

In response to the proposed scenario, the remote MAS succeeded to bring the system back to desired state (figure 7). Besides validating the good functionality and the scalability of the proposed control based on MAS, this is also an interesting demonstration of the application of our proposed consistency condition in combining RTS and HIL from multiple infrastructures in a holistic validation of smart grid.

Conclusion

Co-simulation and Hardware-in-the-loop are two emerging techniques for the assessment and validation of cyber-physical energy system. While there are existing implementations of these technologies in the literature, it is necessary to establish a theoretical framework for the configuration of these experiments. In this paper, we proposed a formalization of the synchronization aspect in co-simulation, hardware-in-the-loop as well as the integration of these two technologies. Based on this formalization, we proposed several conditions in which a co-simulation framework is well-posed and correctly synchronized. We then considered HIL framework as a special case of the synchronization conditions, we established a relaxed realtime constraint for such a set up. On one hand, it allows specifying the parameters of a hardware-in-the-loop framework to satisfy the relaxed RT constraint. On the other hand, it enables the possibility to achieve a reasonable timing consistency with non-real-time simulation, thus, paves the way towards HIL with slow-simulation and variable time steps.

In general, the analysis in this paper provides a theoretical base to evaluate whether a co-simulation and HIL framework is well-posed and consistently synchronized. It emphasizes that the synchronization in a co-simulation framework need to be actualized upon the individual simulation horizons which are not necessarily linear with the micro/marco steps. Same thing goes for HIL experiment as timing consistency is underlined over the speed of the simulator. Remote/multi-platform cosimulation and integration of HIL to co-simulation for more complex CPES assessment should respect and follow the same principles.

The proposed synchronization condition is still strict in the sense that it requires pairwise coherency of the simulation horizons. A more in-depth study on transitivity and composition of the relation will eventually lead to a more simplified/relaxed version of the condition 6, which would be more useful in application. Besides, as S, R and M are not bijective, the mapping from one time horizon to another will be subjected to information loss (i.e. cannot fully recovered a model behavior from received signals). Finding the "loss-less" condition for such mapping is therefore an interesting perspective.

Figure 1 :

 1 Figure 1: Example of mapping of signal in a co-simulation interface

Definition 4 .Definition 5 . 10 )

 4510 Given a co-simulation framework S of n models S 1 , S 2 , . . . , S n , the co-simulation horizon of the framework is defined as The intersection of these simulator's horizons defines the interaction zone:I S = ⋃ i, j=1..n,i≠ j (T S i ∩ T S j ) (Condition2. (Necessary Condition) The co-simulation framework S is well-posed if:

Condition 4 .

 4 (Two-Sided Condition) Two model A and B are synchronized, noted as A B iff

Condition 7 .

 7 (Relaxed Real-time Condition) The simulation S and the hardware-under-test HuT are correctly coupled, iff

Condition 8 .

 8 (Simplified Real-time Condition) The simulation S and the hardware-under-test HuT are correctly coupled, iff

Figure 2 :

 2 Figure 2: Proposed MAS based control architecture.

Figure 3 :

 3 Figure 3: System Under Test.

Figure 4 :

 4 Figure 4: Cross-infrastructure remote HIL test setup.

Figure 5 :

 5 Figure 5: Histogram of latency, fitting to Kernel distribution and its CDF.

Figure 6 :

 6 Figure 6: Histogram of latency, fitting to Gamma distribution and its CDF.

Figure 7 :

 7 Figure 7: Regulation of frequency (a), voltage (b) and power sharing among ESS (s) with the MAS network.

Acknowledgment

This work is supported by the Carnot Institute Energies du Futur under the PPInterop II project and the European Community's Horizon 2020 Programme under project ERIGrid (Grant Agreement No. 654113). The authors would like to express their gratitude to the anonymous reviewers who helped to improve the paper quality with their remarks.


domain: sdv math

text:
F M • Hamelin
email: fhamelin@agrocampus-ouest.fr

Y • Mammeri

Y • Aigu

S E • Strelkov

M A Lewis

Host diversification may split epidemic spread into two successive fronts advancing at different speeds

Host diversification methods such as within-field mixtures (or field mosaics, depending on the spatial scale considered) are promising methods for agroecological plant disease control. We explore disease spread in host mixtures (or field mosaics) composed of two host genotypes (susceptible and resistant). The pathogen population is composed of two genotypes (wild-type and resistance-breaking). We show that for intermediate fractions of resistant hosts, the spatial spread of the disease may be split into two successive fronts. The first front is led by the wild-type pathogen and the disease spreads faster, but at a lower prevalence, than in a resistant pure stand (or landscape). The second front is led by the resistance-breaking type, which spreads slower than in a pure resistant stand (or landscape). The wild-type and the resistancebreaking genotype coexist behind the invasion fronts, resulting in the same prevalence as in a resistant pure stand. This study shows that host diversification methods may have a twofold effect on pathogen spread compared to a resistant pure stand (or landscape): on one hand they accelerate disease spread, and on the other hand they slow down the spread of the resistancebreaking genotype. This work contributes to a better understanding of the

multiple effects underlying the performance of host diversification methods in agroecology.

Introduction

The negative impacts of pesticides on biodiversity and human health require developing more ecological methods to control plant diseases [START_REF] Matthews | Pesticides: health, safety and the environment[END_REF][START_REF] Yang | Enhanced agricultural sustainability through withinspecies diversification[END_REF][START_REF] Wuest | Ecological and evolutionary approaches to improving crop variety mixtures[END_REF]. So far, the main alternative to using pesticides against plant pathogens has been to breed genetically resistant plant varieties or cultivars and deploy them as pure stands [START_REF] Wolfe | The increased use of diversity in cereal cropping requires more descriptive precision[END_REF]. Under these conditions, populations of pathogens often evolve and break down resistance genes after a few years, while a breeding program may require at least a decade [START_REF] Brown | Durable resistance of crops to disease: a Darwinian perspective[END_REF][START_REF] Zhan | Playing on a pathogen's weakness: using evolution to guide sustainable plant disease control strategies[END_REF]. More sustainable control methods will require diversifying genetic resistance over time (cultivar rotations) [START_REF] Bargués-Ribera | Eco-evolutionary agriculture: Host-pathogen dynamics in crop rotations[END_REF][START_REF] Nilusmas | Multiseasonal modelling of plant-nematode interactions reveals efficient plant resistance deployment strategies[END_REF], and/or in space. Spatial diversification of genetic resistance can be achieved either at the field scale (host mixtures) [START_REF] Jeger | Disease progress of non-specialised fungal pathogens in intraspecific mixed stands of cereal cultivars. I. Models[END_REF][START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF][START_REF] Skelsey | Influence of host diversity on development of epidemics: an evaluation and elaboration of mixture theory[END_REF][START_REF] Ohtsuki | Epidemiology and disease-control under gene-for-gene plantpathogen interaction[END_REF][START_REF] Sapoukhina | Spatial deployment of gene-for-gene resistance governs evolution and spread of pathogen populations[END_REF][START_REF] Sapoukhina | Spatially mixed crops to control the stratified dispersal of airborne fungal diseases[END_REF][START_REF] Iacono | Durable resistance to crop pathogens: an epidemiological framework to predict risk under uncertainty[END_REF][START_REF] Fabre | Epidemiological and evolutionary management of plant resistance: optimizing the deployment of cultivar mixtures in time and space in agricultural landscapes[END_REF][START_REF] Mikaberidze | Developing smarter host mixtures to control plant disease[END_REF][START_REF] Rousseau | Virus epidemics, plant-controlled population bottlenecks and the durability of plant resistance[END_REF][START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF] or at the scale of the landscape (field mosaics) [START_REF] Skelsey | Invasion of Phytophthora infestans at the landscape level: how do spatial scale and weather modulate the consequences of spatial heterogeneity in host resistance?[END_REF][START_REF] Fabre | Durable strategies to deploy plant resistance in agricultural landscapes[END_REF][START_REF] Djidjou-Demasse | Mosaics often outperform pyramids: insights from a model comparing strategies for the deployment of plant resistance genes against viruses in agricultural landscapes[END_REF][START_REF] Lof | Achieving durable resistance against plant diseases: scenario analyses with a national-scale spatially explicit model for a wind-dispersed plant pathogen[END_REF][START_REF] Papaïx | Differential impact of landscape-scale strategies for crop cultivar deployment on disease dynamics, resistance durability and long-term evolutionary control[END_REF][START_REF] Rimbaud | Assessing the durability and efficiency of landscape-based strategies to deploy plant resistance to pathogens[END_REF][START_REF] Rimbaud | Mosaics, mixtures, rotations or pyramiding: What is the optimal strategy to deploy major gene resistance?[END_REF][START_REF] Watkinson-Powell | When does spatial diversification usefully maximize the durability of crop disease resistance?[END_REF].

Host mixtures consist in growing several varieties of the same plant species in the same field and at the same time [START_REF] Wolfe | The current status and prospects of multiline cultivars and variety mixtures for disease resistance[END_REF][START_REF] Mundt | Use of multiline cultivars and cultivar mixtures for disease management[END_REF]. Mixtures are or have been used against plant pathogens in various regions of the world [START_REF] Finckh | Cereal variety and species mixtures in practice, with emphasis on disease resistance[END_REF][START_REF] Zhu | Genetic diversity and disease control in rice[END_REF][START_REF] Han | Intercropping of rice varieties increases the efficiency of blast control through reduced disease occurrence and variability[END_REF][START_REF] Reiss | Cultivar mixtures: a meta-analysis of the effect of intraspecific diversity on crop yield[END_REF]. Although host mixtures have long been studied both theoretically [START_REF] Kampmeijer | Epimul, a simulator of foci and epidemics in mixtures, multilines, and mosaics of resistant and susceptible plants[END_REF][START_REF] Jeger | Disease progress of non-specialised fungal pathogens in intraspecific mixed stands of cereal cultivars. I. Models[END_REF] and experimentally [START_REF] Jeger | Disease progress of non-specialised fungal pathogens in intraspecific mixed stands of cereal cultivars. II. Field experiments[END_REF][START_REF] Wolfe | The current status and prospects of multiline cultivars and variety mixtures for disease resistance[END_REF], their design remains to be optimized to be more widely and efficiently used [START_REF] Mikaberidze | Developing smarter host mixtures to control plant disease[END_REF][START_REF] Ben | Improved control of Septoria tritici blotch in durum wheat using cultivar mixtures[END_REF]. So far, most empirical studies of epidemic control through host diversification have taken place at the field scale, but a few exceptions [START_REF] Mundt | Landscape heterogeneity and disease spread: experimental approaches with a plant pathogen[END_REF][START_REF] Papaïx | Influence of cultivated landscape composition on variety resistance: an assessment based on wheat leaf rust epidemics[END_REF] suggest that host diversification at the landscape scale (field mosaics) would also be effective [START_REF] Papaïx | Can epidemic control be achieved by altering landscape connectivity in agricultural systems?[END_REF].

Host diversification is often achieved with resistant and susceptible plants in which resistance is qualitative, meaning that infection either succeeds or fails (as opposed to quantitative resistance, which only partially decreases the success of infection). Qualitative resistance is often conferred by major resistance genes and driven by gene-for-gene interactions [START_REF] Flor | Current status of the gene-for-gene concept[END_REF][START_REF] Milgroom | Population biology of plant pathogens: genetics, ecology, and evolution[END_REF]. Pathogen genotypes can then be classified into two types: the resistance-breaking type, which can successfully infect both resistant and susceptible hosts, and the wild-type, which can successfully infect susceptible hosts only.

Infection sites may be individual plants, or smaller infection units such as the size of a lesion, if the spatial domain considered is a single field [START_REF] Zadoks | Systems analysis and the dynamics of epidemics[END_REF]. At a larger spatial scale, albeit with a lower spatial resolution, infection sites may be regarded as fields if the spatial domain considered is a landscape. Previous authors [START_REF] Mundt | Development of crown rust epidemics in genetically diverse oat populations: effect of genotype unit area[END_REF] introduced the term "genotype unit area" to denote the contiguous area occupied by a single host genotype [START_REF] Papaïx | Can epidemic control be achieved by altering landscape connectivity in agricultural systems?[END_REF]. From now on, we will use the term "infection site" to mean the same thing as "genotype unit area" at the landscape scale. Infection sites correspond to either resistant or susceptible host genotypes. If the spatial domain is a field, resistant and susceptible infection sites form a host mixture. If the spatial domain is a landscape, resistant and susceptible infection sites form a mosaic of pure stands in the landscape.

ACCEPTED MANUSCRIPT / CLEAN COPY

Models can be spatially implicit [START_REF] Jeger | Disease progress of non-specialised fungal pathogens in intraspecific mixed stands of cereal cultivars. I. Models[END_REF][START_REF] Ohtsuki | Epidemiology and disease-control under gene-for-gene plantpathogen interaction[END_REF][START_REF] Fabre | Durable strategies to deploy plant resistance in agricultural landscapes[END_REF][START_REF] Iacono | Durable resistance to crop pathogens: an epidemiological framework to predict risk under uncertainty[END_REF][START_REF] Fabre | Epidemiological and evolutionary management of plant resistance: optimizing the deployment of cultivar mixtures in time and space in agricultural landscapes[END_REF][START_REF] Mikaberidze | Developing smarter host mixtures to control plant disease[END_REF][START_REF] Djidjou-Demasse | Mosaics often outperform pyramids: insights from a model comparing strategies for the deployment of plant resistance genes against viruses in agricultural landscapes[END_REF][START_REF] Rousseau | Virus epidemics, plant-controlled population bottlenecks and the durability of plant resistance[END_REF][START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF] or spatially explicit [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF][START_REF] Skelsey | Influence of host diversity on development of epidemics: an evaluation and elaboration of mixture theory[END_REF][START_REF] Sapoukhina | Spatial deployment of gene-for-gene resistance governs evolution and spread of pathogen populations[END_REF][START_REF] Sapoukhina | Spatially mixed crops to control the stratified dispersal of airborne fungal diseases[END_REF][START_REF] Skelsey | Invasion of Phytophthora infestans at the landscape level: how do spatial scale and weather modulate the consequences of spatial heterogeneity in host resistance?[END_REF][START_REF] Lof | Achieving durable resistance against plant diseases: scenario analyses with a national-scale spatially explicit model for a wind-dispersed plant pathogen[END_REF][START_REF] Rimbaud | Assessing the durability and efficiency of landscape-based strategies to deploy plant resistance to pathogens[END_REF][START_REF] Rimbaud | Mosaics, mixtures, rotations or pyramiding: What is the optimal strategy to deploy major gene resistance?[END_REF][START_REF] Papaïx | Differential impact of landscape-scale strategies for crop cultivar deployment on disease dynamics, resistance durability and long-term evolutionary control[END_REF][START_REF] Watkinson-Powell | When does spatial diversification usefully maximize the durability of crop disease resistance?[END_REF]. See [START_REF] Rimbaud | Models of plant resistance deployment[END_REF] for a recent review. Among spatially explicit models, the few that tackle the speed of pathogen spread focus on host mixtures [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF][START_REF] Skelsey | Influence of host diversity on development of epidemics: an evaluation and elaboration of mixture theory[END_REF]. In these models, the pathogen population is assumed to be wild-type, that is unable to infect resistant hosts. In this study, we are interested in comparing the spreading speeds of resistance-breaking and wild-type pathogens invading a spatial domain where resistant and susceptible infection sites are mixed.

An outline of the paper is as follows. We first present the canola-clubroot pathosystem as a motivating example. Then, we describe the epidemic model in Section 3, along with its biological interpretation. We look into the analogous spatially implicit (diffusionless) system in Section 4, and show how the model parameters drive the genetic composition of the pathogen population at equilibrium. We study the spatial spread of the disease in Section 5: we present formal derivations of the spreading speeds, that are checked with numerical simulations. Finally, in Section 6, we present a discussion of our results.

2 Motivating example: Clubroot in Alberta, Canada Clubroot, caused by the obligate parasite Plasmodiophora brassicae Wor., is a soilborne disease of canola/oilseed rape (Brassica napus L) and other cruciferous hosts. Infection by P. brassicae is associated with the formation of large galls on the roots of susceptible plants, leading to yield losses estimated at 10%-15% globally [START_REF] Dixon | Biology of Plasmodiophora brassicae wor. -A review of recent advances[END_REF]. In Alberta, Canada, clubroot was first identified on canola in 2003 [START_REF] Tewari | Identification of clubroot of crucifers on canola (Brassica napus) in alberta[END_REF] and by 2019 had been detected in more than 3000 fields [START_REF] Strelkov | The occurrence and spread of clubroot on canola in alberta in 2019[END_REF]. Targeted surveys for the occurrence and spread of clubroot have been conducted annually since 2005, generating a large data set. The total infested area is estimated as the number of fields in which clubroot has been confirmed, multiplied by the standard size (160 acres) of an agricultural field in Alberta. Figure 1 shows the spread of clubroot as the square root of infested acres over time (years). The linear relationship indicates an approximately constant speed of disease spread.

The disease is controlled mainly by planting clubroot resistant canola varieties, which initially became available in 2009 and 2010. Resistance-breaking (RB) genotypes of P. brassicae were detected for the first time in Alberta in 2013 [START_REF] Strelkov | Emergence of new virulence phenotypes of plasmodiophora brassicae on canola (brassica napus) in alberta, canada[END_REF]. By 2019, RB genotypes had been confirmed in 239 clubroot-infested fields vs. 2928 fields infested by the wild-type (WT) genotype. The resistancebreaking designation of each P. brassicae isolate is confirmed in a greenhouse test, although given the large number of new cases of clubroot identified every year, only isolates recovered from resistant varieties are tested for their ability to overcome resistance. As above, the spread of the WT and RB genotypes of P. brassicae is expressed, respectively, as the square root of WT and RB infested acres over time (years) (Fig. 2). The spread of the RB genotypes is apparently slower than that of the WT. This suggests the RB genotypes might have a lower capacity to propagate compared to the WT.

Model

We consider a spatial domain where susceptible and resistant infection sites are homogeneously mixed for simplicity.

There are two pathogen genotypes: the wild-type (WT) and the resistancebreaking (RB) genotype. The WT can only infect susceptible hosts, whereas the RB genotype can infect both susceptible and resistant hosts. For the sake of simplicity, we do not allow for coinfections in the model, meaning that each site can be infected by at most one pathogen genotype. This assumption fits resistant infection sites which cannot be infected by the WT. However, this is a strong assumption with respect to susceptible fields at the landscape scale. This is a weaker assumption for susceptible sites at the field scale, depending on whether the infection is systemic (as in viruses) or localized in the plant (as in rust fungi, which form very small lesions that may be defined as infection sites as well).

The total density of infection sites N in the spatial domain is assumed to be homogeneous and constant. A fraction p of the infection sites is resistant and the other fraction 1 -p is susceptible. The density of susceptible and resistant sites infected by the RB genotype are I s and I r , respectively. The density of susceptible sites infected by the WT is J s . Plant pathogens are generally ACCEPTED MANUSCRIPT / CLEAN COPY dispersed as free-living propagules such as fungal spores, or by insect vectors in the case of viruses. For simplicity we do not model possible insect vectors explicitly. The density of infectious propagules produced by the WT and the RB genotype are W and V , respectively. Each site infected by the WT emits infectious propagules at rate ζ.

We assume the RB genotype incurs a cost c which reduces the propagule emission rate ζ by a factor 0 ≤ 1 -c ≤ 1 relative to the WT. Similarly, a cost k reduces the infection efficiency of the propagules. The idea of a cost as a counterpart to the ability of breaking a resistance gene originated as a theoretical hypothesis to explain the often observed maintenance of polymorphism in pathogen populations, in both agricultural and wild ecosystems [START_REF] Vanderplank | Disease resistance in plants[END_REF][START_REF] Sasaki | Host-parasite coevolution in a multilocus gene-for-gene system[END_REF][START_REF] Gandon | The evolution of parasite virulence, superinfection, and host resistance[END_REF][START_REF] Tellier | Stability of genetic polymorphism in host-parasite interactions[END_REF][START_REF] Brown | Durable resistance of crops to disease: a Darwinian perspective[END_REF]. Since then, such a cost has been demonstrated and measured in a number of parasites, including bacteria [START_REF] Cruz | Predicting durability of a disease resistance gene based on an assessment of the fitness loss and epidemiological consequences of avirulence gene mutation[END_REF][START_REF] Wichmann | Effector genes of Xanthamonas axonopodis pv. vesicatoria promote transmission and enhance other fitness traits in the field[END_REF], fungi [START_REF] Carson | Aggressiveness and perennation of isolates of Cochliobolus heterostrophus from North Carolina[END_REF][START_REF] Thrall | Evolution of virulence in a plant host-pathogen metapopulation[END_REF][START_REF] Bahri | Tracking costs of virulence in natural populations of the wheat pathogen, Puccinia striiformis f. sp. tritici[END_REF][START_REF] Huang | Fitness cost of virulence differs between the Avrlm1 and Avrlm4 loci in Leptosphaeria maculans (phoma stem canker of oilseed rape)[END_REF][START_REF] Caffier | Aggressiveness of eight Venturia inaequalis isolates virulent or avirulent to the major resistance gene Rvi6 on a non-Rvi6 apple cultivar[END_REF][START_REF] Bruns | The jack of all trades is master of none: A pathogen's ability to infect a greater number of host genotypes comes at a cost of delayed reproduction[END_REF][START_REF] Bousset | Spatio-temporal connectivity and host resistance influence evolutionary and epidemiological dynamics of the canola pathogen Leptosphaeria maculans[END_REF], viruses [START_REF] Jenner | A fitness cost for Turnip mosaic virus to overcome host resistance[END_REF][START_REF] Janzac | A point mutation in the polymerase of Potato virus Y confers virulence toward the Pvr4 resistance of pepper and a high competitiveness cost in susceptible cultivar[END_REF][START_REF] Fraile | Rapid genetic diversification and high fitness penalties associated with pathogenicity evolution in a plant virus[END_REF][START_REF] Poulicard | Why Rice yellow mottle virus, a rapidly evolving RNA plant virus, is not efficient at breaking rymv1-2 resistance[END_REF][START_REF] Ishibashi | Coevolution and hierarchical interactions of Tomato mosaic virus and the resistance gene Tm-1[END_REF][START_REF] Khatabi | Fitness penalty in susceptible host is associated with virulence of Soybean mosaic virus on Rsv1-genotype soybean: a consequence of perturbation of HC-Pro and not P3[END_REF], nematodes [START_REF] Castagnone-Sereno | Selection and parasite evolution: a reproductive fitness cost associated with virulence in the parthenogenetic nematode Meloidogyne incognita[END_REF], and oomycetes [START_REF] Montarry | Fitness costs associated with unnecessary virulence factors and life history traits: evolutionary insights from the potato late blight pathogen Phytophthora infestans[END_REF].

The diffusion coefficient of infectious propagules is D. Each infectious propagule dies or deposits at rate δ. Their infection rate is θ. Infected sites recover or are replaced with uninfected sites at rate α. We assume no diseaseinduced mortality for simplicity, as is the case for most plant viruses and many other parasites. Time is denoted as t and space is denoted as x ∈ (-∞, +∞) (i.e. we adopt a uni-dimensional conception of space for simplicity).

The model is:

∂W ∂t = ζJ s -δW + D ∂ 2 W ∂x 2 , ∂V ∂t = (1 -c)ζ(I s + I r ) -δV + D ∂ 2 V ∂x 2 , ∂J s ∂t = θW ((1 -p)N -J s -I s ) -αJ s , ∂I s ∂t = (1 -k)θV ((1 -p)N -J s -I s ) -αI s , ∂I r ∂t = (1 -k)θV (pN -I r ) -αI r .

Non-dimensionalization

Let

w = δW ζN , v = δV ζN , a = J s N , y = I s N , z = I r N ,
and

t * = tδ , x * = x δ D ,
and

R = θζN αδ , µ = α δ .
Dropping the asterisks, and using subscripts to denote partial differentiation with respect to t or x, the model can be expressed as:

w t = a -w + w xx , v t = (1 -c)(y + z) -v + v xx , a t = µ[Rw(1 -p -a -y) -a] , (1)
y t = µ[(1 -k)Rv(1 -p -a -y) -y] , z t = µ[(1 -k)Rv(p -z) -z] .

Analysis of the diffusionless system

To evaluate the underlying dynamical behavior, we first explore the diffusionless system, in which the dot represents differentiation over time:

ẇ = a -w , v = (1 -c)(y + z) -v , ȧ = µ[Rw(1 -p -a -y) -a] , (2)
ẏ = µ[(1 -k)Rv(1 -p -a -y) -y] , ż = µ[(1 -k)Rv(p -z) -z] .
ACCEPTED MANUSCRIPT / CLEAN COPY This system of ODE's is similar to that of [START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF], except that we consider explicit propagule compartments, which allows us to distinguish two types of resistance-breaking costs (c and k).

We define P = a + y + z as the prevalence of the disease.

4.1 The system is Cooperative System (2) has "Defined Feedback Relations" [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF] meaning that the rates of change of all variables depends either positively or negatively of other variables (no sign change). More specifically, the rate of change of each variable in (w, v, a, y, z) depends non-negatively on other variables with two exceptions: the rate of change of a depends negatively on y, and the rate of change of y depends negatively on a (sign symmetry). Figure 3 shows the influence graph corresponding to System (2). One can check that every loop in the signed and undirected influence graph has an even number of negative edges. Hence, system (2) is cooperative [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]. Therefore, the solution of the system necessarily converges to an equilibrium (there is no attracting periodic orbit for instance).

Equilibria

System (2) has four equilibria:

the "disease-free" equilibrium (w, v, a, y, z) = (0, 0, 0, 0, 0).

the "WT-only" equilibrium ( w, 0, ā, 0, 0

) := R(1 -p) -1 R , 0, R(1 -p) -1 R , 0, 0 ,
which is positive and therefore biologically exists if and only if

R(1 -p) > 1 . (3)
The prevalence of the disease is

P W = (1 -p) -1/R.

ACCEPTED MANUSCRIPT / CLEAN COPY

the "RB-only" equilibrium (0, v, 0, ȳ, z), in which

v := R(1 -c)(1 -k) -1 R(1 -k) , ȳ := (1 -p) R(1 -c)(1 -k) -1 R(1 -c)(1 -k) , z := p R(1 -c)(1 -k) -1 R(1 -c)(1 -k) .
The RB-only equilibrium is positive and therefore biologically exists if and only if

R(1 -c)(1 -k) > 1 . ( 4
)
The prevalence of the disease is

P R = 1 -1/(R(1 -c)(1 -k)). -the "coexistence" equilibrium ( ŵ, v, â, ŷ, ẑ), in which ŵ := c -p+k(1 -c) c+k(1 -c) , v := Rp(1 -c) -c-k(1 -c)(Rp + 1) Rc+kR((1 -c)(1 -k) + c) , â := c -p+k(1 -c) c+k(1 -c) , ŷ := Rp(1 -c) -c-k(1 -c)(Rp + 1) Rc+kR(1 -c) , ẑ := Rp(1 -c) -c-k(1 -c)(Rp + 1) R(1 -c)(1 -k) .
The coexistence equilibrium is positive and therefore biologically exists if and only if

c+k(1 -c) > p and p > c+k(1 -c) R(1 -c)(1 -k) , (5)
which implies that conditions (3) and ( 4) are satisfied. The prevalence of the disease is again

P R = 1 -1/(R(1 -c)(1 -k)).
Figure 4A summarizes the picture in the plane (p, c) for R = 5 and k = 0. Figure 4B shows the prevalence of the disease (P ) as a function of the fraction of resistant sites (p) for R = 5, c = 0.5, and k = 0. The shape of the graphs are qualitatively unchanged for other parameter values. As previously shown by [START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF], the prevalence decreases with respect to p (the proportion of resistant sites) until a threshold p value is attained, above which the prevalence remains constant.

To make clear under which conditions coexistence occurs, we next perform a mutual invasion analysis [START_REF] Samia | Coexistence and spread of competitors in heterogeneous landscapes[END_REF].

ACCEPTED MANUSCRIPT / CLEAN COPY

Can the WT invade the RB genotype?

First we wonder whether the WT can initially invade the RB-only equilibrium

(assuming R(1 -c)(1 -k) > 1). The sub-model associated with this question is ẇ = a -w , ȧ = µ[Rw(1 -p -a -ȳ) -a] , ACCEPTED MANUSCRIPT / CLEAN COPY with ȳ = (1 -p) R(1 -c)(1 -k) -1 R(1 -c)(1 -k) .
This yields ẇ = a -w ,

ȧ = µ Rw 1 -p R(1 -c)(1 -k)
-a -a .

The corresponding Jacobian matrix is

J (w,a) = -1 1 µR 1-p R(1-c)(1-k) -a -µ(Rw + 1)
, which we evaluate at (w, a) = (0, 0):

J (0,0) = -1 1 µ 1-p (1-c)(1-k) -µ .
The trace of J (0,0) is negative and its determinant is

µ(p -c-k(1 -c)) (1 -c)(1 -k) .
Therefore, the wild-type can invade if and only if

c+k(1 -c) > p . (6)
4.4 Can the RB genotype invade the WT?

Now we wonder whether the RB genotype can initially invade the WT-only equilibrium (assuming R(1 -p) > 1). The sub-model associated with this question is

v = (1 -c)(y + z) -v , ẏ = µ[(1 -k)Rv(1 -p -ā -y) -y] , ż = µ[(1 -k)Rv(p -z) -z] , with ā = R(1 -p) -1 R .
Letting u = y + z, the above system can be equivalently expressed as

v = (1 -c)u -v , u = µ (1 -k)Rv 1 - R(1 -p) -1 R -u -u .

ACCEPTED MANUSCRIPT / CLEAN COPY

The corresponding Jacobian matrix is

J (v,u) = -1 1 -c µ(1 -k)R 1 -R(1-p)-1 R -u -µ((1 -k)Rv + 1) ,
which we evaluate at (v, u) = (0, 0):

J (0,0) = -1 1 -c µ(1 -k)(Rp + 1) -µ ,
The trace of J (0,0) is negative and its determinant is

µ(1 -(1 -k)(Rp + 1)(1 -c)) .
Therefore, the WT can invade if and only if

p > c+k(1 -c) R(1 -c)(1 -k) . (7)
4.5 Conclusions from the diffusionless system

The WT and RB pathogen can invade the disease-free equilibrium if and only if conditions (3) and (4) hold, respectively, which we assume (otherwise, we would not be concerned with these pathogen genotypes).

The WT can initially invade the RB-only equilibrium if condition (6) holds. The RB genotype can initially invade the WT-only equilibrium if condition (7) holds. Therefore, if both conditions ( 6) and ( 7) hold, the RB and WT can invade each other, and they are led to coexist. Conditions ( 6) and ( 7) are the same as conditions (5) implying the existence of the coexistence equilibrium. Since system (2) is cooperative and other equilibria are unstable, the coexistence equilibrium is globally asymptotically stable when it exists. We refer the reader to [13, Section S4] for more details on how the theory of cooperative systems can be used to show the global stability of the coexistence equilibrium in a related model.

Spreading of interacting pathogen genotypes with multiple fronts

We now get back to the spatial model (with diffusion). Although we are not able to rigorously prove the existence of travelling wave fronts, our numerical evidence strongly suggests that such fronts exist (Figure 5). Therefore we consider the possibility of travelling fronts for our system. Figure 5 shows the six possible patterns of travelling waves. We posit that, as can be found in some similar cooperative degenerate systems [START_REF] Li | Traveling wave solutions in partially degenerate cooperative reaction-diffusion systems[END_REF][START_REF] Fang | Traveling waves for monotone semiflows with weak compactness[END_REF] (see also [START_REF] Hamelin | Mate finding, sexual spore production, and the spread of fungal plant parasites[END_REF][START_REF] Doli | Phénomènes de propagation de champignons parasites de plantes par couplage de diffusion spatiale et de reproduction sexuée[END_REF]), the front speeds are linearly determined as given as minimum possible wave speed based on the linearization at the leading edge of the wave. In this section, we apply the minimum wave speed approach [START_REF] Lewis | Biological invasion of an organism with separate mobile and stationary states: Modeling and analysis[END_REF][START_REF] Hadeler | Spatial dynamics of the diffusive logistic equation with a sedentary compartment[END_REF][START_REF] Bampfylde | Biological control through intraguild predation: case studies in pest control, invasive species and range expansion[END_REF][START_REF] Janson | Resultant and discriminant of polynomials[END_REF] to linearized submodels for finding the pathogen spreading speeds as a critical point in each sub-model, and then test the theoretical results against numerical simulations.

ACCEPTED MANUSCRIPT / CLEAN COPY

Invading the disease-free equilibrium

In this section, we derive the spreading speeds of the WT and RB genotype when invading the disease-free equilibrium (Fig. 5A-B).

At leading edge invading the disease-free equilibrium, w, v, a, y, z have small positive values. Note that w and v equations in system (1) are linear. We ACCEPTED MANUSCRIPT / CLEAN COPY linearize the other equations at leading edge:

w t = a -w + w xx , v t = (1 -c)(y + z) -v + v xx , a t ≈ µ[Rw(1 -p) -a] , y t ≈ µ[(1 -k)Rv(1 -p) -y] , z t ≈ µ[(1 -k)Rvp -z] .
Using u = y + z, the above system can be equivalently expressed as

w t = a -w + w xx , a t = µ[Rw(1 -p) -a] , v t = (1 -c)u -v + v xx , u t = µ[(1 -k)Rv -u] ,
where (w, a) and (v, u) now form two uncoupled subsystems.

We are interested in traveling wave solutions such that

=     w a v u     = k exp(-s(x -σt)) ,
in which k is an implicit column vector, σ is the wave speed, and s is the exponential decay rate of the wave profile at leading edge.

Plugging the previous expression in the system, we get

sσ =     -1 + s 2 1 0 0 µR(1 -p) -µ 0 0 0 0 -1 + s 2 (1 -c) 0 0 µ(1 -k)R -µ     A , which implies det     -1 + s 2 -sσ 1 0 0 µR(1 -p) -µ -sσ 0 0 0 0 -1 + s 2 -sσ (1 -c) 0 0 µ(1 -k)R -µ -sσ     B=A-sσ1 = 0 ,
in which 1 is the identity matrix. The above equality is equivalent to

det -1 + s 2 -sσ 1 µR(1 -p) -µ -sσ × det -1 + s 2 -sσ (1 -c) µ(1 -k)R -µ -sσ = 0 , ACCEPTED MANUSCRIPT / CLEAN COPY which is equivalent to 0 = [(-1 + s 2 -sσ)(-µ -sσ) -µR(1 -p)] × [(-1 + s 2 -sσ)(-µ -sσ) -µR(1 -c)(1 -k)] .
We now introduce Q to denote either R(1 -p) or R(1 -c)(1 -k), which are the basic reproductive numbers of the WT and RB genotypes, respectively; see equations [START_REF] Bampfylde | Biological control through intraguild predation: case studies in pest control, invasive species and range expansion[END_REF][START_REF] Bargués-Ribera | Eco-evolutionary agriculture: Host-pathogen dynamics in crop rotations[END_REF]. From now on, we therefore assume Q > 1 so that the pathogens can invade the disease-free equilibrium. This allows us to consider a single equality:

(-1 + s 2 -sσ)(-µ -sσ) -µQ = 0 . ( 8
)
Next we follow the approach for using this equation to calculate the minimum wave speed as outlined in [START_REF] Hadeler | Spatial dynamics of the diffusive logistic equation with a sedentary compartment[END_REF]. Solving equation ( 8) for σ > 0 yields

σ(s, µ, Q) := s 2 -µ -1 + (s 2 + µ -1) 2 + 4µQ 2s .
Equation ( 8) can therefore be written to include the dependency of σ on s, µ and Q as

P(σ(s, µ, Q), s) := (-1 + s 2 -sσ(s, µ, Q))(-µ -sσ(s, µ, Q)) -µQ = 0 .
Differentiating with respect to s, we have, for all s,

dP ds = ∂P ∂σ ∂σ ∂s + ∂P ∂s = 0 . (9)
We are interested in the minimum possible wave speed. Let

s (µ, Q) = arg min s σ(s, µ, Q) , and
σ (µ, Q) = σ(s (µ, Q), µ, Q) .
Since σ is such that ∂σ/∂s = 0, equation ( 9) yields

∂P ∂s (σ (µ, Q), s (µ, Q)) = 0 . ( 10
)
Since P is cubic in s, ∂P/∂s is quadratic in s. We are interested in the conditions on the coefficients that allow both polynomials to have a common root, s . They are given by cancelling the resultant of the two polynomials (which is also the determinant of the Sylvester matrix of the two polynomials). Letting P = es 3 + f s 2 + gs + h yields ∂P/∂s = 3es 2 + 2f s + g, and the resultant is

r(e, f, g, h) = -e(f 2 g 2 -4eg 3 -4f 3 h + 18ef gh -27e 2 h 2 ) ,
as described in equation (4.30) of reference [START_REF] Janson | Resultant and discriminant of polynomials[END_REF]. The coefficients are identified as e = -σ, f = σ 2 -µ, g = σ(1 + µ), and h = -µ(Q -1). The equality r(e, f, g, h) = 0 can be equivalently expressed as a cubic with respect to σ 2 : in which

E(σ 2 ) 3 + F (σ 2 ) 2 + G(σ 2 ) + H = 0 , ACCEPTED MANUSCRIPT / CLEAN COPY
E = 4µ(Q -1) + (1 + µ) 2 , F = 4 + 2µ 3 + (6Q + 2)µ 2 + (18Q -8)µ, G = µ 4 -(6Q -8)µ 3 -(27Q 2 -36Q + 8)µ 2
, and H = -4µ 4 (Q -1). Since we assume Q > 1, we have that E is positive and H is negative, which means that we are in the same configuration as [START_REF] Hadeler | Spatial dynamics of the diffusive logistic equation with a sedentary compartment[END_REF]. This implies that σ (µ, Q) is uniquely defined as the square root of the largest root of the above cubic.

Although it is possible to write down the formula for the largest root of a cubic polynomial, we have no simple expression of σ (µ, Q). However, this function only depends on two parameters, µ and Q. Therefore, one can plot σ as a function of µ and Q without loss of generality. One can easily check that for all Q > 1 and µ > 0,

∂σ ∂Q (µ, Q) > 0 , (11)
meaning that σ increases as Q increases. As a consequence, c+k(1 -c) > p , or equivalently (1

-p) > (1 -c)(1 -k), implies σ 1 := σ (µ, R(1 -p)) > σ (µ, R(1 -c)(1 -k)) =: σ 2 ,
meaning that the WT spreads faster than the RB genotype. For c+k(1 -c) < p, the inequality is reversed, and the RB genotype spreads faster than the WT.

Invading the RB-only equilibrium

In this section, we derive the spreading speed of the WT when invading the RB-only equilibrium (Fig. 5C-E). At leading edge invading the RB-only equilibrium, w, a have small positive values while v, y, z are constant (equal to equilibrium values v, ȳ, z, resp.). We ACCEPTED MANUSCRIPT / CLEAN COPY linearize the a equation at leading edge:

w t = a -w + w xx , a t ≈ µ[Rw(1 -p -ȳ) -a] , with ȳ = (1 -p) R(1 -c)(1 -k) -1 R(1 -c)(1 -k) .
This yields

w t = a -w + w xx , a t ≈ µ 1 -p (1 -c)(1 -k) w -a .
Proceeding as in Section 5.1, the spreading speed of the WT is

σ 3 := σ µ, 1 -p (1 -c)(1 -k) .
As shown in Section 5.1, σ 3 > 0 is conditioned to the second argument of the function σ being greater than one, i.e. (1

-p)/[(1 -c)(1 -k)] > 1.
This condition is equivalent to [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF], which is the condition for the WT to invade the RB genotype. Lastly, one can check that σ 3 < σ 1 is equivalent to condition (4), which is necessarily satisfied since we assumed the RB-only equilibrium exists.

Invading the WT-only equilibrium

In this section, we derive the spreading speed of the RB genotype when invading the WT-only equilibrium (Fig. 5D-F). At leading edge invading the WT-only equilibrium, v, y, z have small positive values while w, a are constant (equal to equilibrium values w, ā, resp.) in space. We linearize the (y, z) equations at leading edge:

v t = (1 -c)(y + z) -v + v xx , y t ≈ µ[(1 -k)Rv(1 -p -ā) -y] , z t ≈ µ[(1 -k)Rvp -z] , with ā = R(1 -p) -1 R .
Using u = y + z, the above system can be equivalently expressed as

v t = (1 -c)u -v + v xx , u t = µ[(1 -k)(1 + Rp)v -u] .
ACCEPTED MANUSCRIPT / CLEAN COPY Proceeding as in Section 5.1, the spreading speed of the RB genotype is

σ 4 := σ (µ, (1 + Rp)(1 -c)(1 -k)) .
As indicated in Section 5.1, σ 4 > 0 is conditioned to the second argument of the σ function being greater than one, i.e. (1

+ Rp)(1 -c)(1 -k) > 1.
This condition is equivalent to [START_REF] Bousset | Spatio-temporal connectivity and host resistance influence evolutionary and epidemiological dynamics of the canola pathogen Leptosphaeria maculans[END_REF], which is the condition for the RB genotype to invade the WT. One can lastly check that σ 4 < σ 2 is equivalent to condition (3), which is necessarily satisfied since we assumed the WT-only equilibrium exists. In the remaining situation, that is

c+k(1 -c) > p > c+k(1 -c) R(1 -c)(1 -k) ,
the WT invades the disease-free equilibrium faster than the RB genotype (since c+k(1 -c) > p) but does not exclude the later behind the invasion front. Therefore, the RB genotype invades the WT-only equilibrium at speed σ 4 < σ 1 (since σ 4 < σ 2 and σ 2 < σ 1 in this case). Behind this second and slower invasion front, the WT and the RB genotype coexist steadily (Fig. 7). The reverse situation, in which the RB genotype would precede the WT while both are led to coexist, is not possible.

ACCEPTED MANUSCRIPT / CLEAN COPY

Note that the spreading speed of the RB genotype invading the disease-free equilibrium (σ 2 ) does not depend on p, so if the WT is faster than RB in a mixture, it is also faster than RB in a pure stand.

Numerical simulations

Numerical computations of system (1), presented in Figure 5 and Figure 7, are performed using the fourth order Runge-Kutta method for the first derivative with respect to time, combined with the second order finite difference scheme for the second derivative with respect to space.

Figure 7 shows the long-term pattern obtained from numerical simulations with parameters such that both pathogen genotypes (WT and RB) can coexist.

If the WT initially precedes the RB genotype, the WT spreads faster than the RB genotype and generates a first invasion front. For these parameter values, the RB genotype can still invade and coexist with the WT. This leads to the emergence of a second invasion front, advancing more slowly than the first front. Some authors refer to stacked fronts, e.g. [START_REF] Iida | Stacked fronts for cooperative systems with equal diffusion coefficients[END_REF], or propagating terraces, e.g. [START_REF] Polácik | Planar propagating terraces and the asymptotic one-dimensional symmetry of solutions of semilinear parabolic equations[END_REF]. Between the first and the second front, the prevalence of the disease reaches a plateau. The prevalence of the disease is increased behind the second front, and is eventually the same as in a resistant pure stand (or landscape), see Fig. 4 and supplementary videos.

If the RB genotype initially precedes the WT, the WT makes up the delay and overtakes the RB genotype. In the long-run, the WT invasion front precedes the RB invasion front as in Fig. 7.

Discussion

Developing new or preexisting methods based on bio-diversification forms the basics of agroecology [START_REF] Altieri | Agroecology: the science of sustainable agriculture[END_REF]. Plant host diversification (cultivar mixtures or field mosaics) are one of these. Several models explored the spreading speed of plant diseases in host mixtures [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF][START_REF] Skelsey | Influence of host diversity on development of epidemics: an evaluation and elaboration of mixture theory[END_REF]. In particular, it was shown that the spreading speed is an increasing concave function of the proportion of susceptible hosts [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF], which is consistent with the spreading speed we obtained for the WT invading the disease-free equilibrium. However, previous studies did not consider pathogen diversity, specifically the possibility that the WT co-occur with a RB genotype. Since pathogen diversity and the competition between pathogen genotypes for susceptible hosts is a key component of the success of host mixtures [START_REF] Ohtsuki | Epidemiology and disease-control under gene-for-gene plantpathogen interaction[END_REF][START_REF] Mikaberidze | Developing smarter host mixtures to control plant disease[END_REF][START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF], we explored a simple spatial model accounting for pathogen diversity.

The model we explored consists of 2 PDE's coupled to 3 ODE's. We showed that the model is cooperative, which led us to make a formal spreading speed analysis around each possible equilibrium. This allowed us to get a rather complete picture of the behaviour of the model. In particular, we were able to derive the spreading speed of the RB genotype invading an area previously invaded by the WT.

ACCEPTED MANUSCRIPT / CLEAN COPY

A notable feature of our analysis is that in host mixtures, the disease may spread along two fronts, with maximum prevalence behind the second front. The first front is led by the WT invading the disease-free equilibrium (as in [START_REF] Van Den Bosch | Focus expansion in plant disease. IV: Expansion rates in mixtures of resistant and susceptible hosts[END_REF]), and the second front (advancing at a lower speed) is led by the RB genotype invading the WT. Both pathogen genotypes eventually coexist at the rear of the front. The first front advances faster than the disease spreading speed in a resistant pure stand (or landscape). This is because the resistance-breaking cost exceeds the proportion of resistant hosts in this case. However, the prevalence just behind the first invasion front is lower than in a pure resistant stand (or landscape). This is because the WT is unable to infect resistant hosts. The second invasion front advances more slowly than the disease spreading speed in a resistant pure stand (or landscape). This is because the spread of the RB genotype is slowed down by the WT due to the competition for susceptible hosts. Behind the second invasion front however, the prevalence is the same as in a resistant pure stand (or landscape). Altogether, our results show that the spatial spread of the disease in host mixtures (or field mosaics) may have a twofold effect compared to its spread in a resistant pure stand (or landscape): on one hand the WT invades the disease-free area faster than the RB would, and on the other hand the spread of the RB genotype is slowed down due to the competition with the WT.

Our study therefore showed that the question of whether host mixtures (or field mosaics) are advantageous compared to resistant pure stands (or landscapes) is not straightforward, be it in terms of spreading speed or prevalence. However, we did not consider other features of host mixtures that are known to be key for their performance in the field, such as priming-induced crossprotection [START_REF] Clin | Taking advantage of pathogen diversity and immune priming to minimize disease prevalence in host mixtures: a model[END_REF]. Moreover, we ignored co-infections of susceptible plants (or fields) by the WT and the RB genotype [START_REF] Watkinson-Powell | When does spatial diversification usefully maximize the durability of crop disease resistance?[END_REF]. This model was a first step in these directions, that are left for future research.

Fig. 1

 1 Fig. 1 Spread of clubroot (Plasmodiophora brassicae) in Alberta, Canada, expressed as the square root of P. brassicae-infested acres over time. The total number of infested acres was calculated by multiplying the number of fields with confirmed clubroot infestations by 160 (standard acreage of an agricultural field in Alberta). The linear regression equation is shown on the figure.

Fig. 2

 2 Fig. 2 Spread of wild-type (WT) and resistance-breaking (RB) genotypes of the clubroot pathogen Plasmodiophora brassicae in Alberta, Canada, expressed as the square root of WT (black) and RB (red) infested acres over time. The total number of infested acres was calculated by multiplying the number of fields with confirmed infestations of each genotype by 160 (standard acreage of an agricultural field in Alberta). The linear regression equations are shown in black and red, respectively, for the WT and RB data.

Fig. 3

 3 Fig. 3 Influence graph corresponding to System (1).

Fig. 4

 4 Fig. 4 Parameter values: R = 5, c = 0.5, and k = 0.

Fig. 5

 5 Fig. 5 Six possible patterns of travelling waves. Panel (A) shows a travelling wave solution connecting the disease-free equilibrium (DFE) to the WT-only equilibrium. Panel (B) shows the connection betwenn the DFE and the RB-only equilibrium. Panel (C) shows the connection between the RB-only and the WT-only equilibrium. Panel (D) shows the symmetric situation in which the RB genotype invades the WT. Panel (E) shows the connection between the RB-only and the Coexistence equilibrium. Panel (F) shows the connection between the WT-only and the Coexistence equilibrium. The invasion fronts spread to the right at speeds σ 1 (A), σ 2 (B), σ 3 (C, E), and σ 4 (D, F). Parameter values: p = 0.2, c = 0.9 (A), p = 0.9, c = 0.2 (B), p = 0.25, c = 0.7 (C), p = 0.5, c = 0.25 (D), p = 0.3, c = 0.5 (E), p = 0.2, c = 0.25 (F), and R = 5, µ = 0.1, k = 0 (A-F).

Fig. 6

 6 Fig.6Value of σ with respect to Q and µ, which is non positive if and only Q ≤ 1. Note that σ has no biological meaning for Q ≤ 1.

Fig. 7

 7 Fig. 7 Numerical simulations showing multiple fronts. See also the supplementary videos. Parameter values: R = 5, c = 0.5, p = 0.4, µ = 0.1, k = 0. Initial conditions: w(x, 0) = v(x, 0) = a(x, 0) = 0.5, y(x, 0) = 0.25, z(x, 0) = 0.25 for 0 ≤ x ≤ 20, and 0 elsewhere..

5. 4

 4 Comparing spreading speeds under coexistence conditions If c+k(1 -c) < p, then the RB genotype invades the disease-free equilibrium faster than the WT and competitively excludes the later behind the invasion front. If c+k(1 -c) > p and p < [c+k(1 -c)]/[R(1 -c)(1 -k)], then the reverse situation happens: the WT invades the disease-free equilibrium faster than the RB genotype and competitively excludes the later behind the invasion front.

ACCEPTED MANUSCRIPT / CLEAN COPY

Acknowledgements FH acknowledges funding from the INRAE "Plant Health and the Environment" Division. MAL gratefully acknowledges a Canada Research Chair and an NSERC Discovery grant. The authors thank the reviewers for their helpful suggestions.

Data Availability Statement

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest

The authors declare that they have no conflict of interest.