Modelling temperature-dependent dynamics of single and mixed infections in a plant virus

Applied Mathematical Modelling 102 (2022) 694–705 Contents lists available at ScienceDirect Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm Modelling temperature-dependent dynamics of single and mixed infections in a plant virus Josep Sardanyés a , b , ∗, Cristina Alcaide c , Pedro Gómez c , Santiago F. Elena d , e a Centre de Recerca Matemàtica (CRM). Edifici C, Campus de Bellaterra, 08193 Cerdanyola del Vallès, Barcelona, Spain b Dynamical Systems and Computational Virology, CSIC Associated Unit CRM-Institute for Integrative Systems Biology (I 2 SysBio), Spain c Centro de Edafología y Biología Aplicada del Segura (CEBAS), CSIC, Departamento de Biología del Estrés y Patología Vegetal, PO Box 164, 30100, Murcia, Spain d I 2 SysBio, CSIC-Universitat de València, Parc Científic UV, Paterna 46182 València, Spain e Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA a r t i c l e i n f o Article history: Received 4 June 2021 Revised 1 October 2021 Accepted 4 October 2021 Available online 13 October 2021 Keywords: Abiotic stress Bifurcations Co-infection dynamics Dynamical systems Nonlinear dynamics Thermal reaction norms Transcritical bifurcations a b s t r a c t Multiple viral infection is an important issue in health and agriculture with strong impacts on society and the economy. Several investigations have dealt with the population dynamics of viruses with different dynamic properties, focusing on strain competition during multiple infections and the effects on viruses’ hosts. Recent interest has been on how multiple infections respond to abiotic factors such as temperature ( T ). This is especially important in the case of plant pathogens, whose dynamics could be affected significantly by global warming. However, few mathematical models incorporate the effect of T on parasite fitness, especially in mixed infections. Here, we investigate simple mathematical models incorporating thermal reaction norms (TRNs), which allow for quantitative analysis. A logistic model is considered for single infections, which is extended to a Lotka-Volterra competition model for mixed infections. The dynamics of these two models are investigated, focusing on the roles of T -dependent replication and competitive interactions in both transient and asymptotic dynamics. We determine the scenarios of co-existence and competitive exclusion, which are separated by a transcritical bifurcation. To illustrate the applicability of these models, we ran singleand mixed-infection experiments in plants growing at 20 ◦C and 30 ◦C using two strains of the plant RNA virus Pepino mosaic virus . Using a macroevolutionary algorithm, we fitted the models to the data by estimating the TRNs for both strains in single infections. Then, we used these TRNs to feed the mixedinfection model estimating the strength of competition. We found an asymmetrical pattern in which each strain dominated at different T values due to differences in their TRNs. We also identified that T can modify competition interference greatly for both isolates. The models proposed here can be useful for investigating the outcomes of multiple-infection dynamics under abiotic changes and have implications for the understanding of viral responses to global warming. ©2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) ∗Corresponding author. E-mail address: jsardan[email protected] (J. Sardanyés). https://doi.org/10.1016/j.apm.2021.10.008 0307-904X/© 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 1. Introduction Many examples exist of different parasites infecting a host simultaneously [1,2] . This is of special importance because multiple infection can cause huge impacts on health and agriculture, thus having severe ecological and socio-economic consequences [3,4] . Regarding diseases impacting human health, human immunodeficiency virus type-1 (HIV-1) [5,6] can coinfect with tuberculosis (TB) [7] , the hepatitis B [8,9] and C viruses [10] , and malaria [11] . Other common examples of multiple infection include infection by the hepatitis B and C viruses [12] , gonorrhea and chlamydia [13] , and herpes simplex viruses 1 and 2 [14,15] . In all cases, the outcomes of multiple infection differ from the observed from the single infection cases. The dynamics of different pathogenic microbial strains infecting the same host mainly have been studied with dynamical systems theory such as for general diseases [16–19] and microparasites (viruses, bacteria, protozoa, or fungi) [20–28] . More recently, several studies have focused on the dynamical outcomes of different evolving virus strains with different infectious phenotypes, i.e. , specialist vs. generalist, infecting host cells [29,30] . Multiple infection dynamics, in the form of both coinfections and superinfections, have also been an object of intense mathematical modeling [22,23] . Mixed viral infection is especially relevant for plant crops [2,31–34] . In this sense, within-plant virus–virus interactions can affect a plant’s epidemiology as a result of synergistic or antagonistic interactions among infecting strains [33] or viral species [2] . They also have major consequences on virulence and virus fitness [34–37] . Emerging plant viruses are responsible for widespread crops edpidemics, representing a major challenge to plant health and thus to agriculture [38–40] . This emergence has driven by intrinsic viral and host factors, in addition to ecological, agronomic and socio-economic factors [41–44] . An important yet poorly explored issue is how environmental abiotic stresses such as droughtiness, salinity, light intensity, and temperature alter plant physiology and thus affect coinfecting viruses at the within-host scale. This becomes an important question given the ongoing global warming scenario [45] . Climate change is likely to increase the frequency of viral diseases emerging in plant crops [46–48] . Warming and highly variable climate may directly and indirectly affect host, vectors, and viral traits, further influencing viral epidemics in both wild and cultivated plants. Therefore, exploring how and to what extent temperature variations may affect the eco-evolutionary dynamics of viral populations can contribute to a more complete understanding of how global warming will affect virus emergence and epidemiology. Typically, multi-strain mathematical models use different fitness traits, e.g. , replication capacity, to determine how competitive interactions affect transients and equilibrium dynamics [25–30,49] . However, as mentioned above, the successful replication and infection of these pathogens may also depend on environmental abiotic factors. That is, different strains may have different responses to such environmental components of the ecosystem. In this article, we introduce and investigate mathematical dynamical models explicitly considering temperature ( T )- dependent replication. We do so by considering the so-called thermal reaction norm (TRN). Briefly, reaction norms relate fitness with some environmental variable: the TRN describes how a particular virus strain responds to different temperatures. The model for single infections considers a time-continuous logistic model, which describes exponential replication (at small population values) and includes intra-specific competition, together with the T -dependent replication. This simple model is extended to a Lotka–Volterra competition model with two strains. The dynamics of these models are investigated analytically and numerically, focusing on the role of T -dependent replication. The two viruses competition model reveals a transition between coexistence and out-competition through a transcritical bifurcation. As a practical application, we finally use these models to fit experimental data gathered for Pepino mosaic virus (PepMV; genus Potexvirus , family Alphaflexiviridae ) strains infecting tomato plants at 20 ◦C and 30 ◦C, considering both single and mixed infections (see also [50] ). PepMV has a major impact on tomato production worldwide. PepMV populations mainly comprise a mixture of two types of co-circulating isolates: the European (EU) and Chilean (CH2) strains [51] . The CH2 type isolates are predominant. However, EU isolates have not been displaced but persist mainly in mixed infections. The CH2 strain has higher fitness than the EU strain does [51] ; hence, the question remaining to be answered is what mechanisms may promote coexistence, despite this fitness advantage of the CH2 strain. Our model for mixed infection confirms the stable coexistence of both strains under different fitness values, while also allowing both the TRNs and competition coefficients to be quantified. 2. Mathematical models In the next sections, we introduce and analyze the mathematical models for single and mixed infection at different temperatures ( T ). The models assume well-mixed viral populations and are aimed at describing within-host ( in planta ) virus replication. As mentioned at the end of the Introduction section, we used a logistic time-continuous model for single infection, including a T -dependent replication (fitness) parameter. We will refer to the viral load with the variable x i . Because the experimental results were obtained using the EU and CH2 strains, we will hereafter use subindices Eand Cto represent them, respectively. 2.1. Modelling the effect of temperature in virus replication In this section, we introduce the function used to model T -dependent replication rates, r i (T ) with i ∈ { E, C} , i.e. , the TRN. This function (see [52–54] for details) reads as follows: r i (T ) = a i (T −T min )  1 −exp (b i (T −T max ))  . (1) 695 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Fig. 1. Thermal reaction norms modeled by Eq. (1) , which provide replication rates, r i (T ) , as a function of temperature. We show results for a i = 0 . 2 (a), a i = 0 . 4 (b), and a i = 0 . 6 (c), using different values of b i : 0.010 (black), 0.025 (red), 0.050 (green), 0.075 (blue), and 0.100 (orange). Panel (d) was computed by setting a i = 0 . 6 and using larger values of b i : 0.10 (black), 0.50 (red), 0.75 (green), 1.00 (blue), and 1.50 (orange). Equation (1) is a bell-shaped function assuming that the PepMV genome’s (gRNA) replication response to T is nonlinear, increasing as T grows until reaching an optimal T for replication. These nonlinear thermal traits have improved the predictions for several diseases (see [45] for a review). Parameters a i and b i allow the bell’s height and shape to be tuned. As our first approach, a i mainly determined the magnitude of the optimal T while b i defined the asymmetry of the TRN, thus influencing the location of the optimal T . In our analyses, we set the range of temperatures between T min = 15 ◦C ≤T ≤50 ◦C = T max . To illustrate the behavior of Eq. (1) , we show, in Fig. 1 , the replication values for different values of a i and b i . As an application, constants a i and b i will be estimated in Section 3.2.1 from the experimental data for single infection to determine the replication rates of each strain at the two studied temperatures, thus obtaining the TRN for each PepMV strain. 2.2. Dynamical equations for single and mixed infections The dynamics of virus replication for single infection was modeled using the following one-variable model: f(x i ) = dx i dt = r i (T ) x i 1 −x i K . (2) Here, variable x i > 0 is the viral load of each strain i ∈ { E, C} . This equation is the well-known logistic model, here including a T -dependent replication rate (TRN function given by Eq. (1) ). The growth rate is exponential for small population sizes, and the virus populations will grow until the carrying capacity, K > 0 (ng gRNA), is achieved. The dynamics for mixed infection are modeled similarly to single infection is modeled, while also considering T - dependent replication rates and including competition between both strains. The strength of the competition by strain j on strain i is introduced with (dimensionless) parameters βij > 0 , i, j ∈ { E, C} , i  = j. We use the classic Lotka–Volterra competition model, given by the following: dx E dt = r E (T ) x E 1 −x E + βEC x C K , (3) dx C dt = r C (T ) x C 1 −βCE x E + x C K , (4) where the two viruses compete for the same cellular resources, Kbeing the carrying capacity. 696 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Fig. 2. (a) Dynamics (shown in linear-log scale) obtained from Eq. (5) for 60 days post-inoculation (dpi) with K = 10 4 , a i = 0 . 2 , and both b i = 0 . 01 (black) and b i = 0 . 1 (orange). As the initial conditions, we used x i (0) = 1 . The solid and dashed lines correspond to temperatures of 20 ◦C and 30 ◦C, respectively. The inset displays the phase space of Eq. (2) , with two equilibrium points: P ∗ 1 (white circle) and P ∗ 2 (blue circle). The arrows indicate the stability of both equilibria. (b) The potential function obtained from Eq. (6) using the same parameter values as used in panel (a) and represented with the same colors and line styles. The local minimum corresponds to the equilibrium point P ∗ 2 = K, which is stable. 3. Results and discussion In this section, we first discuss the dynamics of the models for both single and mixed infection, focusing on the role of T -dependent replication. Then, we estimate the parameters describing the TRN of each PepMV strain from single-infection experiments. Finally, we quantify the competition coefficients for both strains from the experimental mixed infection. An overview of the experimental protocols used to characterize the single and double infections with the CH2 and EU strains is provided in Sections A1 and A2 of the Appendix. 3.1. Dynamics Equation (2) can be solved analytically, giving x i (t) = x i (0) K x i (0) + (K −x i (0)) exp (−r i (T ) t) , (5) where x i (0) > 0 is the initial conditions and tis time. Despite being a well-known model, let us summarize the dynamical properties of the logistic model, now considering the T -dependent replication rate. This system has two equilibrium points, namely P ∗ 1 = 0 , and P ∗ 2 = K. The stability of this equilibrium is determined from λ= df(x i ) dx i = r i (T ) 1 −2 x i K . It is easy to show that P ∗ 1 is unstable and P ∗ 2 is an attractor for T min < T < T max since λ(P ∗ 1 ) = r i (T ) > 0 and λ(P ∗ 2 ) = −λ(P ∗ 1 ) . That is, for any initial condition, the virus population will always achieve the carrying capacity. Here, the equilibrium values do not depend on r i (T ) but on how they are approached or left do so. We want to emphasize that the carrying capacity of plant cells might vary depending on the temperature conditions; thus, different equilibria may be achieved at different temperatures (as observed in the experiments). Hence, the values of Kwill be estimated from the experiments (see Section B below). The time dynamics of Eq. (2) can be visualized as shown in Fig. 2 (a), by means of time series obtained from the explicit solution given by Eq. (5) . The inset in this figure displays a schematic diagram of the phase space, with the equilibria P ∗ 1 being unstable and P ∗ 2 being stable. As mentioned, T has an important effect on the duration of the transients since the eigenvalues are hyperbolic and proportional to r i (T ) . For the parameter values chosen in Fig. 2 (a), x i achieves equilibrium faster at 30 ◦C (dashed lines) that at 20 ◦C (solid lines). Another way to visualize the dynamics in one-variable dynamical systems is by computing the potential function, here given by U(x i ) = − f(x i ) dx i = r i (T ) x 2 i x i 3 K −1 2 . (6) The potential is a cubic polynomial function and appears parabolic in the studied range of x i ( Fig. 2 (b)), with a minimum at the equilibrium value P ∗ 2 = K = 10 4 . Notice that the wells are deeper at 30 ◦C, meaning that the equilibrium point is much more attracting, in turn explaining why this equilibrium is achieved faster than at 20 ◦C. However, the dynamics will be slower for T s above the peak in the TRN (see next section). 697 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 We sough to examine competition by studying the dynamics of Eqs. (3) - (4) , which have four different equilibrium points, namely P ∗ 0 = (0 , 0) , P ∗ E = (K, 0) , P ∗ C = (0 , K) , and P ∗ EC = (x ∗ E , x ∗ C ) , with x ∗ E = K(βEC −1) (βEC βCE −1) , x ∗ C = K(βCE −1) (βEC βCE −1) . Note that the fixed points are also independent of the replication rates, r i (T ) , now also being determined by Kand the competition constants. However, as stated below, the replication rates also play an important role in the transients toward equilibria. The Jacobian matrix reads J = r E (T )(K −2 x E −x C βEC ) /K −r E (T ) x E βEC /K −r C (T ) x C βCE /K r C (T )(K −2 x C −x E βCE ) /K . The eigenvalues are here given by λ1 (P ∗ 0 ) = r E (T ) , λ2 (P ∗ 0 ) = r C (T ) , and thus P ∗ 0 is a repeller since r i (T ) > 0 . The eigenvalues of equilibrium point P ∗ E are λ1 (P ∗ E ) = −r E (T ) , λ2 (P ∗ E ) = r C (T )(1 −βCE ) . Notice that the first eigenvalue is negative (when r E (T ) > 0) , and the stability of this equilibrium depends on the second eigenvalue. For r E (T ) > 0 , r C (T ) > 0 , and βCE > 1 , the equilibrium P ∗ E is an attractor. For r E (T ) > 0 , r C (T ) > 0 and βCE < 1 the equilibrium P ∗ E is a saddle point. The stability of the equilibrium P ∗ C is similar to the explained above. Here, the eigenvalues are: λ1 (P ∗ C ) = r E (T )(1 −βEC ) , λ2 (P ∗ C ) = −r C (T ) . This equilibrium will be stable when βEC > 1 and r E (T ) > 0 and will be a saddle when βEC < 1 and r E (T ) > 0 (since the second eigenvalue is always negative). The eigenvalues for equilibrium P ∗ EC are λ±(P ∗ EC ) = r E (T )(1 −βEC ) + r C (T )(1 −βCE ) ∓√  2 ( βEC βCE −1) , with = (r E (T ) +r C (T )) 2 +4 r E (T ) r C (T ) (βEC βCE −1) , = βEC −1 and = βCE −1 . A global picture of the dynamics is provided in Fig. 3 . Panel (a) displays the outcompetition of CH2 by EU, while panel (b) displays the strains’ coexistence (as we observed in the experiments performed in this work). Panel (c) provides a bifurcation diagram plotting the equilibrium values of both strains at increasing values of βEC . Notice that for βEC < 1 , the two viral strains coexist. At βEC > 1 , the coexistence turns into EU being outcompeted by CH2 (see the phase portrait in the inset). The panel below the bifurcation diagram shows the eigenvalues for the fixed points P ∗ C and P ∗ EC . Both fixed points collide and interchange stability by means of one eigenvalue through a transcritical bifurcation. Finally, Fig. 4 shows the temporal dynamics, obtained numerically (numerical simulations were performed with a fourth-order Runge–Kutta method with a constant time step t = 0 . 05 .) for different temperature values. Here, as shown for the single-infection dynamics, temperatures involving higher replication rates of viral gRNA will allow the equilibrium points to be achieved earlier. However, for temperatures close to the maximum values ( e.g. , 40 ◦C and 49 ◦C) the dynamics slow down again. 3.2. Application of the models to the experimental data In this section, we use the previously described mathematical models to estimate parameters from the experimental data for single and mixed infection with EU and CH2 isolates. Appendix Sections A1 and A2 provide information about the experiments. Further details about experiments can be found in Ref. [50] , in which the evolutionary outcomes for both single and mixed infection were studied using next-generation sequencing. The experimental data processing for the models’ fitting is commented on in Section A3 The optimization method used for the fittings, given by a macroevolutionary algorithm (MA), is described in Section A4 The MA is used to compute the parameters providing the best fits after running 50 replicates, as well as the mean values of the optimized parameters and of the best fittings averaged over these replicates. 3.2.1. Single infections For the single-infection experiment, we used Eq. (2) to estimate the model parameters, considering the vector (a i , b i , K) , i ∈ { E, C} , which includes the parameters of the TRN and the plants’ carrying capacities. Here, we show the values of the parameters providing the best fit from all 50 replicates of the MA. The mean values ( ±SD ) of the optimized parameters and of the best parameter sets averaged over the replicates are displayed in supplementary Table SI, together with the best parameter sets, as commented on below. 698 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Fig. 3. Dynamics of Eqs. (3) - (4) . (a) Phase portrait showing the outcompetition of x C by x E , with βEC = 0 . 5 and βCE = 1 . 5 (the arrows indicate the direction of the orbits; the equilibrium points are shown with solid circles: repeller [white], saddle point [grey], attractor [black]). (b) Coexistence dynamics with βEC = 0 . 25 . Coexistence was observed in the experiments discussed in Section 3.2.2. (c) Bifurcation diagram increasing the competition strength by x C on x E , computed with 0 . 05 ≤βEC ≤1 . 95 and βCE = 0 . 5 . Notice that the two strains coexist for 0 < βEC < 1 . The solid line and the solid dots in black indicate the equilibrium values for x E obtained numerically and analytically, respectively. In red, we show the same results for x C . At βEC = 1 , equilibria P ∗ EC and P ∗ C collide in a transcritical bifurcation, and for βEC > 1 , the equilibrium P ∗ C becomes the attractor. The lower panel displays the eigenvalues for the fixed points P ∗ C and P ∗ EC . In all of the analyses, a E = a C = 0 . 8 , b E = b C = 0 . 02 , K = 10 4 , and T = 20 ◦C. Fig. 4. Time dynamics (days post-inoculation [dpi] in linear-log scale) for the viral strains EU ( x E (t) in panel (a)) and CH2 ( x C (t) in panel (b)) at different temperatures: 17 ◦C (black); 20 ◦C (red); 25 ◦C (blue); 30 ◦C (green); 35 ◦C (orange); 40 ◦C (violet); and 49 ◦C (grey). Here, we have used a E = 0 . 2 , b E = 0 . 1 , βCE = 0 . 65 , a C = 0 . 4 , b C = 0 . 05 , βEC = 0 . 5 , and K = 10 4 . The parameters providing the best fit for EU at 20 ◦C were a E = 1 . 5519 , b E = 0 . 8138 , r E = 7 . 7595 ng gRNA/day, and K = 1 . 9682 ×10 6 ng gRNA (black trajectory in Fig. 5 a), gRNA being genomic RNA. The parameters providing the best fit for the EU data at 30 ◦C were a E = 1 . 2681 , b E = 0 . 3841 , r E = 19 . 0123 ng gRNA/day, and K = 4 . 115 ×10 5 ng gRNA (red trajectory in Fig. 5 a). Concerning CH2, the parameters providing the best fit at 20 ◦C were a C = 1 . 1382 , b C = 0 . 0151 , r C = 2 . 0757 ng gRNA / day , and K = 4 . 110 ×10 6 ng gRNA . This fitting is displayed in Fig. 5 (b) with a black line. The parameters giving the best fit for 30 ◦C experiments were given by a C = 3 . 5948 , b C = 1 . 3388 , r C = 53 . 9226 ng gRNA/day, and K = 5 . 159 ×10 6 ng gRNA (the fitting using these parameters is displayed in Fig. 5 (b) in red). The fittings of the same ex699 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Fig. 5. Fitting of the mathematical model to the experimental data using the best vector of parameters obtained from 50 replicates of the macroevolutionary algorithm (MA) for single infection of the EU (a) and CH2 (b) strains. The open circles show experimental data at 20 ◦C (black) and 30 ◦C (red); the time trajectories show the dynamics obtained with the mathematical model; and the solid dots show the initial conditions. The fittings were performed using the explicit solution given by Eq. (5) . perimental data above, but obtained from the mean values of the optimized parameters and of the best parameter fits, are shown in supplementary Fig. S1 using the same colors as in Fig. 5 . The replication rates for the EU strain appeared to be slightly higher at 30 ◦C, although the equilibrium value was lower than at 20 ◦C, indicating a slower initial growth of EU at 20 ◦C, due to its carrying capacity being larger at this T . Regarding CH2 dynamics, its carrying capacities were larger than those of EU, with replication rates faster at 30 ◦C than at 20 ◦C, as described for EU. Figure S2 displays how optimization was performed using the MA. Specifically, it shows the evolution of the mean distances,  d i  −, along the generations of the MA, computed from the first, i.e. , the best N/ 4 vectors of parameters during the optimization for the EU (panel (a) at 20 ◦C; panel (b) at 30 ◦C) and CH2 (panel (c) at 20 ◦C; panel (d) at 30 ◦C) strains (Fig. S2). The main panels show the overlapping results for five replicates. The insets display the decrease of the distances for the best vectors of parameters at each improvement event, which are shown overlapped for 25 replicates. The estimation of parameters a i and b i from single-infection experiments allows the reconstruction of the TRN for each strain. According to the experimental data and Eq. (1) , the replication rates for both strains follow a similar TRN. The increase of replication rates with T follows a linear behavior until T gets close to the maximum value, at which point it experiences a sharp decline. These values should be taken with caution since virus replication at high temperatures could vary due to physiological changes in plant cells as well as different enzymatic activities of molecular components that are crucial for gRNA amplification and virus assembly. At intermediate temperatures, the replication values may likely follow this linear fashion. Fig. 6 displays these results for strains EU (a) and CH2 (b). Each plot shows six predictions: three for the experiments at 20 ◦C (solid lines) and three for the ones performed at 30 ◦C (dashed lines). Here the green curves display the profiles obtained from the a i and b i values for the best vector of parameters obtained from the full 50 replicas of the MA. The black curves have been computed with the mean values of a i and b i obtained at the end of the MA algorithm and averaged over the 50 replicates. Finally, the blue lines show the profiles for the values of a i and b i averaged over the best vectors obtained for each replicate of the MA. Regardless the values of a i and b i (especially for the best vector of CH2 at 30 ◦C), the shapes of the estimated TRNs remain similar, meaning that this result, despite the limitations of the experimental data and the selected function, r i (T ) , remains consistent. Future research may experimentally check whether virus replication falls into this type of TRN and locate the optimal T value above which virus replication slows down as T increases. 3.2.2. Mixed infections In the previous section, we estimated the parameter values for the single infections, obtaining the T -dependent replication rates and the carrying capacities for each viral strain at each studied temperature. The main goal of the experiments with the mixed infections was to characterize the competition between EU and CH2 strains at different tem peratures. To do so, we fixed the values of r E (T ) and r C (T ) using the constants a E , b E , and a C , b C which gave the best fits in the single700 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Fig. 6. Viral thermal reaction norms (TRNs) for each isolate, obtained from the estimated parameters of the experimental data of Section 3.2.1. Curves r E (T ) for EU (a) and r C (T ) for CH2 were built from the best parameters (green); from the mean values of the optimized parameters (black) for the 50 replicas of the MA; and from the mean of the best vectors (blue), also computed from the replicas of the MA. The solid and dashed lines show the data for the experiments at 20 ◦C and 30 ◦C, respectively. Fig. 7. Fitting of the experimental data (open circles) for mixed infections (EU strain [black], CH2 strain [green]) obtained with the mathematical model (solid time series) given by Eqs. (3) - (4) using the parameters obtained with the best fit from 50 replicas of the MA. Panels (a) and (b) display the dynamics of competition at T = 20 ◦C and T = 30 ◦C, respectively. The same results for the mean values of the optimized parameters and for the best parameter vectors averaged over the replicates are shown in Fig. S3. infection experiments. Because the two viral strains were co-infecting the same cells in this experiment, the carrying capacity values estimated in the single-infection might differ. Hence, we will introduce the vector (βEC , βCE , K) as free parameters, with βEC and βCE being the competition coefficients among the virus strains. Notably, the competition terms introduced in the mathematical model affect the net replication rate of the viral strains. As mentioned, in order to obtain clear results about the interference due to competition, we used the estimated replication rate of each viral strain at each temperature: r E (T = 20) = 7 . 7594 , r C (T = 20) = 2 . 0757 , r E (T = 30) = 19 . 0122 , and r C (T = 30) = 53 . 9225 . The parameter values providing the best fits for the competition experiments performed at 20 ◦C were βEC = 0 . 9319 , βCE = 0 . 0793 , and K = 7 . 0418 ×10 6 . The results for the experiments at 30 ◦C were βEC = 0 . 9291 , βCE = 0 . 749 , and K = 1 . 4967 ×10 6 . The mean values obtained for the optimized and best parameter vectors averaged over the 50 replicates are also displayed in Table 1 . The time dynamics obtained from the parameters providing the best fits are shown in Fig. 7 . While the CH2 strain exerted a stronger competitive effect at the two studied temperatures, the EU had stronger competition at 30 ◦C than at 20 ◦C. These results suggest that temperature can greately modify the competitive interference between these two isolates, resulting in important changes at the levels of within-host transient dynamics and equilibrium values. 701 J. Sardanyés, C. Alcaide, P. Gómez et al. Applied Mathematical Modelling 102 (2022) 694–705 Table 1 Parameter values estimated using the MA for the competition experiments at 20 ◦C and 30 ◦C using Eqs. (3) - (4) . The best parameters displayed in the first three rows of the table were used to fit the experimental data from Fig. 7 . Parameters Competition at 20 ◦C Competition at 30 ◦C best βEC 0 . 9319 0.9291 best βCE 0 . 0793 0.7497 best K7 . 0418 ×10 6 ng gRNA 1 . 4967 ×10 6 ng gRNA ¯ βEC ±SD 0 . 9319 ±1 . 8 ×10 −7 0 . 9290 ±4 . 16 ×10 −5 ¯ βCE ±SD 0 . 0793 ±5 . 9 ×10 −8 0 . 7490 ±7 . 37 ×10 −4 ¯ K ±SD 7 . 0418 ×10 6 ±13 . 28 ng gRNA 1 . 4 94 9 ×10 6 ±982 . 47 ng gRNA  best βEC  ±SD 0 . 9318 ±1 . 4 ×10 −6 0 . 9291 ±4 . 17 ×10 −5  best βCE  ±SD 0 . 0793 ±4 . 9 ×10 −7 0 . 7490 ±7 . 38 ×10 −4  best K  ±SD 7 . 0418 ×10 6 ±60 . 92 ng gRNA 1 . 4 94 8 ×10 6 ±1050 . 09 ng gRNA 4. Conclusions In this article, we combined mathematical models incorporating temperature-dependent replication rates for viral strains in single and mixed infections with experimental research using two strains of Pepino mosaic virus (PepMV): the European (EU) and the Chilean (CH2) strains. The single-infection model allowed us to estimate the relationship between the viral replication and the temperature of growth of the host i.e. , the thermal reaction norm of the two strains . Thermal reaction norms of parasites are scarce in the literature, which has largely focused on host species [55] . By means of the twodimensional Lotka–Volterra model, we explored the dynamics of strain’s competition, together with the T -dependent replication. We provided a linear stability analysis of the equilibrium points and identified, as expected, a transcritical bifurcation separating the coexistence phase with out-competition. Then, we used the fitness parameters estimated using data with the single-infection model to quantify the strength of interference at increasing temperatures during mixed infections. Interestingly, the coexistence equilibria were independent of replication rates, instead being determined by the carrying capacity and interference coefficients. However, the linear stability analysis indicated that the (temperature-dependent) replication rates indeed affected the transients. We found that the CH2 strain interfered more strongly at 20 ◦C. Our models, despite their simplicity, may be useful for future studies relating intrinsic ecological dynamics (such as competition) to changes in T , and may be of interest to model mixed virus dynamics under future climatic scenarios. Data accessibility . The raw experimental data can be obtained from Dryad at: https://datadryad.org/stash/share/ xTPRVHo8V1LSFWBJFPZ2USy3pipcsggkCIlLg2vAy6s . Authors’ contributions . PG conceived of and designed the study; JS and SFE conceived of the mathematical model; JS analyzed the mathematical models and programmed the optimization algorithm; CA performed the experiments; and all of the authors participated in data analysis, wrote the article, and gave final approval for publication. Funding . J.S. has been partially funded by the CERCA Program of the Generalitat de Catalunya , by the Ministry of Economy, Industry and Competitiveness (MINECO) grant MTM2015-71509-C2-1-R, by Agencia Estatal de Investigación (AEI) grant RTI2018-098322-B-I00, and by a Ramón y Cajal contract (RYC-2017-22243). C.A. was funded by the MINECO within a PhD program grant (FPU16/02569). This work was also supported by the AEI-FEDER grants AGL2014-59556-R and AGL201789550-R to P.G. and PID2019-103998GB-I00 to S.F.E. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. Declaration of Competing Interest We declare we have no competing interests. Appendix A A1. Plant growth conditions, virus inoculation, and sampling The viral stocks used in all of the experiments were produced as follows. Two infectious PepMV clones belonging to the EU [56] and the CH2 [51,57] strains were used to agroinfiltrate Nicotiana benthamiana Domin plants. At 14 days postinoculation (dpi), viral particles were purified from the homogenized plant tissue, following a series of centrifugations and precipitation [58,59] . Tomato plants ( Solanum lycopersicum L. cv. Money Maker) were grown in a greenhouse with a photoperiod of 16 h light: 8 h dark and a temperature between 22 ◦C and 26 ◦C. At 30 days post-germination, two sets of 33 plants were placed into a greenhouse with temperature conditions of 20 ◦C or 30 ◦C. For each temperature, six tomato plants were mock-inoculated, and nine plants per treatment were mechanically inoculated with the EU strain, the CH2 strain, or an equimolar mixture of both strains, i.e. , mixed infections. The inoculations were performed on the third and fourth true leaves by rubbing carborundum and a suspension of virions particles at 500 ng/ μL in sodium phosphate buffer (30 mM). 702