Introduction

Viruses transport their genetic material inside a protein shell, the capsid, surrounded in some species by a lipid membrane. In most viral species, each viral particle contains all the genetic material needed to carry out replication inside a host cell, and generate a progeny of viral particles. A prominent exception to this behavior is found in multipartite viruses. These viruses, first described in the 1960s, have a genome segmented in two or more parts. According to current evidence, the segments are encapsidated separately and, apparently, propagate independently [2, 3]. As of today, there is no mainstream theory able to explain the adaptive advantage of such a strange lifestyle. The main puzzle regarding multipartite viruses is how the simultaneous presence of multiple segments, which imposes severe constraints on the number of viral particles that have to reach a susceptible host, is balanced by other adaptive advantages of multipartitism, whether microscopic or ecological [5, 6].

Despite this apparent paradox, multipartitism is widespread in the Virosphere, as up to 40% of all known viral families are multipartite. A large majority of them infects either plants or fungi, with only four known examples of species infecting exclusively animals. Evolutionary pathways leading to multipartitism are likely to be multiple, since this strategy is present in RNA and DNA viruses, and in the latter case an origin to a single ancestral virus cannot be traced. Beyond its virological interest, multipartitism has a particularly negative effect on agricultural production, as several multipartite viruses are pathogenic, and routinely cripple crop yield. Cultivars themselves may have directly played a role in the rise of multipartitism, as an evolutionary radiation in the diversity of viral species, many being just centuries old, was likely promoted by an intensification of agricultural practices [9–11]. It has been put forward that multipartite species might be at an advantage in the face of environmental changes, since they likely adapt faster due to new combinations of segments promoted by their genomic architecture. It is known that changes in land cover offer multiple opportunities for novel interactions between plants and pathogens [12–15]. Studies on the impact of agriculture in viral ecology have uncovered a surprisingly negative association between plant diversity and family-level diversity of plant-associated viruses, and a higher prevalence of viruses in cultivated areas.

The emergence of defective variants out of the wild-type (wt) form (i.e., the one containing the full genome) has been both posited and observed in controlled environments, arising from replication errors and thriving under conditions that ensure high multiplicity of infection (MOI) [18, 19]. Specifically, it has been shown in vitro that two defective forms spontaneously generated by foot-and-mouth disease virus (FMDV has an unsegmented genome formed by ssRNA of positive polarity) can complement each other and quickly substitute the wt form. This strategy has been formally explored in models of competitive dynamics between the wt and a number of complementing segments, which implemented different advantages that could compensate for the cost of an increased MOI. The model mimicked the experimental setting where, in particular, host-cell availability corresponded to that of a well-mixed system. Two of the advantages implemented had been theoretically proposed in the past, though, as of yet, have not received empirical support (a faster replicative ability and a slower accumulation of deleterious mutations in shorter segments [23, 24]) while, in the case of FMDV, it was shown that capsids containing shorter genomes enjoyed a larger average lifetime between infection events. This differential degradation, dependent on genome length, was sufficient to compensate for co-infection requirements in multipartite forms with two, to up to four, segments, but cannot explain the emergence of multipartite viruses with many segments, such as nanoviruses or babuviruses. Hence, the evolutionary pathway explored in that work would be applicable to a subset of all currently described multipartite viral species.

What is missing from this picture is investigating how the interaction between viral dynamics and host ecology shapes the rise and persistence of multipartite viral forms at the host population level. We also wish to quantify the impact of different host contact structures in driving the success of multipartitism. We tackle this problem by building a compartmental model for studying the competition between monopartite and multipartite variants in terms of their ability to spread and persist on a structured host population.

As the generation of functional defective mutants from the wt occurs at a much longer time scale than the spread of the virus in the population, we set up a model that already contains both the wt and a cohort of defective forms, potentially complementary. This allows us to study the competition dynamics between the different forms causing them to coexist in, or take over the ecological niche. Using both analytical calculations and numerical tools, we investigate the outcome of a random emergence of mutants, and derive the conditions that make multipartitism a fitness-enhancing strategy, allowing the virus to adapt to a wider range of hosts and environments. Since no apparent structural feature discriminates multipartite viruses from monopartite ones—they are found exhibiting different capsid structures, genome sizes and types—we include in the model as few virological features as possible, and investigate how multipartitism impacts on the spreading potential of the virus. We do, however, account for key viral mechanisms that can drive the resilience of multipartitism in an ecological context. The first one is the already mentioned differential degradation, i.e., the different average lifetime of defective viral particles with respect to wt’s, or formally equivalent advantages of faster replication or elimination of deleterious mutations through sex. A second biological mechanism is the mode of transmission of multipartite viruses between hosts. Most known multipartite viruses are spread by vectors (mainly insects), which typically pick up very few viral particles from an infected plant. The transmission process between hosts typically acts as a population bottleneck for the virus, entailing a loss of genetic diversity and, if severe enough, the systematic purge of deleterious forms [28, 29]. Thanks to our parsimonious modeling setup, any of the aforementioned mechanisms can be seen as effectively impacting the chances of the wt or defective particles to reach the target hosts, leading to a difference in transmissibility. A single model parameter, therefore, by tuning this relative transmissibility, embraces a number of different biological processes. In this sense, the results of our model can be extended to other systems as long as the specific mechanisms involved in their spreading fitness can be cast in the form of changes in transmissibility. Remarkably, we find out that even in the absence of an explicit microscopic advantage, ecological dynamics might cause the fixation of the multipartite form due to the stochastic extinction (analogous to random drift) of the monopartite virus.

Alongside the biological properties of the virus, the model implements the structure of contacts among susceptible hosts through which the viruses can spread. Often, in our context, this means the contact network induced by vector movements among plants. Its topology may be diverse, depending on plant distribution and vector behavior, with two limit cases being the distribution of plant species in the wild (see, e.g. and references therein) and huge modern agriculturally homogeneous regions. These different architectures are implemented by tuning the distribution of contact rates among hosts, with limiting cases being fixed contact rate, and power law-distributed contact rate. This feature is the key tool to uncover how the structure of contacts among hosts shapes the endemicity and prevalence of the different viral forms.

We remark that compartmental models of interacting diseases have been studied in the past [32–38]. Those models, however, assume that the disease agents involved are fully-fledged pathogens that can spread on their own, and cannot describe asymmetric viral associations, of which multipartitism is an example. Thanks to that, a completely new phenomenology emerges, driven by a complex evolutionary dynamics, and involving a wide range of ecological interactions: competition, symbiosis, commensality.

The model

We consider a large population of N hosts susceptible to a virus that may circulate in its wild-type (wt) form together with v defective variants that are potentially complementing, i.e., when simultaneously present in the host, they are able to complete the virus infective cycle even in the absence of the wt. As previously done, we take advantage of the timescale separation between the random emergence of defective mutants due to errors during replication of the complete genome, and the competition between forms driven by their spread in the host population. This allows us to effectively model the random emergence of mutants as an initial small, yet nonzero, prevalence in the host population.

As customary in compartmental models, we consider that hosts are either free of the virus, and thus susceptible (S), or infected by a certain combination of the viral forms, translating into various infectious compartments (Fig 1). The main assumption is that a host can be infected only by a combination that guarantees the presence of the full genome. Without it, there is no completion of the viral cycle, and thus no systemic infection is possible. Moreover, we assume that host cells replicate all, and only, the viral forms they are infected by (Fig 1A). These assumptions determine the set of existing compartments.

Two infectious compartments are present regardless of the value of v. They are ⌜wt⌟ and ⌜all⌟, and correspond to plants infected by the wt only, and by the wt together with all the v variants, respectively. If v = 1, no other compartments exist. If v > 1, ⌜seg⌟ identifies plants infected by all the v defective and complementing variants, without the wt. In addition, there are 2v − 1 other compartments containing wt plus a combination of some (not all) of the defective variants. We name them according to which of the latter they are infected by. For instance, ⌜1⌟ contains wt plus variant 1, and ⌜3,5⌟ contains wt plus variants 3 and 5. If defective variants are not present, ⌜wt⌟ behaves like a standard Susceptible–Infected–Susceptible (SIS) model, with probability of transmission upon contact equal to λ.

We implement enhanced transmissibility of defective variants by assuming they spread with a probability ρλ. ρ = 1 thus means that wt and segments are epidemiologically equivalent, while any value larger than 1 causes the defective variants to transmit more easily than the wt. For a graphic representation of the spreading routes and differential transmissibility see Fig 1D. In the general case, an agent in a given compartment may transmit some (or all) of the viral species it hosts to the one it is in contact with, with a probability depending on the initial compartments of the two agents, and on the final compartment. For instance, a host in compartment ⌜1⌟, upon contact with a susceptible one, may transmit both wt and 1, turning the susceptible into a ⌜1⌟. It may instead transmit wt alone, turning the susceptible into a ⌜wt⌟. The defective variant, however, cannot be transmitted alone, as it requires wt, as previously stated. A schematic representation of the compartmental model for a bipartite virus (v = 2) is depicted in Fig 1C.

Our assumption of constant host population (of size N) holds for strictly constant size, as well as populations that are at equilibrium, i.e., the number of births equals the number of deaths, or at least any growth pattern occurs at time scales much larger than the spread of the virus. This assumption is connected to the spreading model, as the recovery process of the Susceptible-Infectious-Susceptible model can be regarded in two ways. It can be seen as proper recovery, with the host clearing the virus but acquiring no immunity to reinfection. It can also be interpreted as the virus killing the host, and a new (susceptible) host filling its ecological space.

In the absence of specific evidence, we make the simplest assumption for transmission: the different types of viral particles are transmitted independently, so that the probability of concurrent transmission of two variants (and wt) is simply the product of the probabilities of the single events. We also assume that co-infection by wt and variants does not alter the infectious period, allowing us to model recovery at a rate μ for all infected hosts. In Methods and S1 Text we expand our analysis to account for nonindependent transmission, heterogeneous recovery rates, and the case when different variants compete for a limited carrying capacity within the host, due, for instance, to a limited number of viral particles a host cell can make per unit time. We show that all these additional features do not impact the qualitative behavior of our model, in agreement with what was previously found in.

When the virus is introduced into a susceptible population, it can either die out quickly and leave the system disease-free (disease-free state, dfs), or reach endemicity. There are four possible endemic states, depending of which variants circulate. We equivalently use the term equilibria, as they are the stable equilibrium points of the spreading dynamics. The first one is wt, in which only the wt is prevalent, and the defective variants have died out. This case maps into an effective SIS model for the compartment ⌜wt⌟. In the second one, hj, any defective segment can circulate alongside the wt because, roughly speaking, the transmissibility of the latter is so high that any defective variant can hijack it, with no need to complement the genome with other variants. In this case, we will likely see the circulation of a number of variants lower than v, as segments can go extinct without hampering the circulation of the remaining ones. The third endemic state, seg, witnesses the presence of all the v segmented variants without wt, and in this case complementation is essential. This state is an SIS model for the compartment ⌜seg⌟. Finally, the state all exhibits circulation of the wt plus all the variants v. Borrowing some terminology from physics, we can then define different epidemic phases. Phases are regions of the space of model parameters. Inside each phase, the macroscopic behavior of viral spread is qualitatively the same. Specifically, we can define a phase in terms of which endemic states it allows. The parameter surfaces separating different phases are called phase transitions, the most important in epidemiology being the epidemic threshold. Below the epidemic threshold, only the dfs exists. Above it, the pathogen can circulate. There, we identify five other phases: wt-phase allows only wt; contingent-phase allows wt, hj and all; mix-phase allows wt, seg and all; seg-phase allows only seg; finally all-phase admits all the possible endemic states.

Table 1 provides a schematic representation of the relationship between phases and endemic states.

In the following, we analytically derive the critical surfaces that separate the different phases in the space of the parameters. This means that, given specific values of the parameters, the possible outcome of the spread can be predicted, thus characterizing the conditions leading to the persistence of multipartitism, and its nature. Then, using numerical simulations, we study the equilibrium prevalences of the endemic states, and their probability of occurring, for a representative set of parameter values.

Firstly, however, we need to set up the theoretical modeling framework in terms of reaction-diffusion equations. For a generic v, we order the compartments by increasing number of viral species they contain, starting from ⌜wt⌟, and ending with ⌜seg⌟, ⌜all⌟. For instance, for v = 3, this would be ⌜wt⌟, ⌜1⌟, ⌜2⌟, ⌜3⌟, ⌜12⌟, ⌜13⌟, ⌜23⌟, ⌜seg⌟, ⌜all⌟. Within the framework of heterogeneous mean field [41–43], we divide hosts in classes according to their contact potential (degree in the language of networks), so that if two hosts have degree k, h, respectively, their contact rate will be the product kh (in the absence of degree-degree correlations). We assume hosts with the same degree are equivalent, and consider the prevalence per degree class. To this end, we define the variable xνk as the prevalence of compartment with index ν and degree class k, i.e., the fraction of the host population which has that degree, and finds itself in that compartment. In terms of xνk, the equations describing the evolution of the disease are

x˙νk=−μxνk+k〈k〉∑β[Γνβ(1−∑σxσk)+∑σΛνβσxσk](∑hhpγ(h)xβh),(1)

where pγ(k) is the probability of a host having degree equal to k. We consider the homogeneous case, where all hosts have the same degree, so that pγ(k) = δk,1 (with no loss of generality we set it to 1), and a highly heterogeneous case, where pγ(k) = Cγ

k−γ is a power-law with exponent γ, and normalization constant Cγ. We denote 〈km〉 as the m-th moment of the degree distribution, computed as 〈km〉 = ∑k

pγ(k)km, as usual. The term 〈k〉 appearing in Eq 1 is then the expected degree. The Greek indices β, ν, σ, run on all the infectious compartments defined before. The susceptible compartment is not included, as the number of susceptible hosts is completely determined by the other compartments, thanks to the assumption of constant population size. Γνβ is the rate of the transition ⌜β⌟⌜S⌟ → β⌜ν⌟, i.e., a transition affecting the prevalence of compartment ⌜ν⌟ through a contact between a host in compartment ⌜β⌟ and a Susceptible. Λνβσ encodes transmission rates among infected individuals, and specifically a transmission from ⌜β⌟ to⌜σ⌟, that leads to the change of the prevalence of ⌜ν⌟. The entries of Γνβ, Λνβσ are functions of λ, ρ and v. Eq (1) thus links the change in the number of hosts with a given degree, and in a compartment (x˙νk), to one reaction and two diffusion processes. The first term, μxνk, represents the decrease due to hosts recovering back to the susceptible state. The second one, with coupling constant Γ, contains the probability of a host, with degree k, being susceptible (1-∑σxσk), and being infected by a host in compartment β, and degree h. The last term, with coupling constant Λ, has the same structure, but the target compartment is a generic infectious compartment σ, instead of the susceptible one. Both infection terms contain the term khpγ(h)/〈k〉, which is the probability a host of degree k establishes a contact with a host of degree h, given a network with no degree-degree correlations.

The analytical approach to computing the critical surfaces consists in studying the linear stability of the different equilibria of Eq (1). Instead, in order to compute the prevalence values and occurrence probability of these equilibria, we have to resort to stochastic spreading simulations. The extensive calculations are reported in Methods and S1 Text, as well as the explanation of the numerical simulations.

Critical behavior

The five phases are completely determined by three surfaces with tractable analytical expressions. They are T1, above which the wt can spread on its own (epidemic threshold for the compartment ⌜wt⌟ while alone); T2, above which segments circulate by hijacking the wt, and Ts, which is the epidemic threshold for the compartment ⌜seg⌟ circulating alone. The expressions we find are

T1={λ=μ^};(2)

T2={λ=1+ρρμ^1+μ^};(3)

Ts={λ=μ^1/vρ}.(4)

μ^ is an effective recovery-rate embodying both the actual recovery rate, and the topology of the contacts: μ^=μ〈k〉/〈k2〉. This entails an important scaling: recovery rate and topology never impact the critical points on their own, but always jointly as μ^. This fact was well-known in the case of the epidemic threshold, Eq (2). Here, we rigorously prove that it extends also to all other critical points. Given that homogeneous contact networks have 〈k〉 ∼ 〈k2〉, the heterogeneous (power-law-like) network recovers the homogeneous case when γ → ∞. Hence, the smaller the exponent γ, the more heterogeneous the contact network is, i.e., hosts with a large number of connections become more likely. These hosts can reach a significant part of the population, and when infectious, they act as superspreades. They are responsible for causing μ^ to go to zero (μ^→0) in the limit of large population size (N → ∞), when the exponent of the degree distribution respects 2 < γ < 3. This implies not only that T1 goes to zero, as it is well-known, but that T2 and Ts do it as well. However, while T2/T1 remains finite as μ^ goes to 0 (T1 and T2 go to zero at the same speed), we find that Ts/T1 → ∞. This entails that Ts goes to zero more slowly, and increasingly so for higher v.

The study of Eqs (2)–(4) reveals four regimes. For low or zero differential transmissibility (ρ<μ^-(v-1)/v-(1-μ^1/v)), as λ increases, one crosses the wt-phase, then the contingent-phase and finally the all-phase (see Fig 2A). For intermediate values of differential transmissibility (μ^-(v-1)/v-(1-μ^1/v)<ρ<μ^-(v-1)/v), the mix-phase substitutes the contingent-phase. This can be seen in Fig 2B and 2C. Then, when μ^-(v-1)/v<ρ<μ^-1, increasing λ causes the system to be in the seg-phase, followed by the mix-phase and later by the all-phase (see Fig 2B, 2C and 2D). Finally, for very high differential transmissibility ρ>μ^-1, the wt no longer spreads and the only possible phase is the seg-phase (see Fig 2B and 2C).

Endemic prevalences

For any possible value of the parameters, Eqs (2)–(4) tell us which endemic states are possible, i.e., which prevalences are higher than zero. They provide, however, no information about the values of such prevalences, which are, in principle, the solutions of the algebraic system obtained by setting x˙ν=0 in Eq (1). A closed-form solution of this system does not exist for heterogeneous networks. In the homogeneous case, while a complete analytical derivation of the endemic states is not possible, we can obtain two important results. Firstly, we notice that the total prevalence of the wt, i.e., the fraction of hosts infected by it (z = ∑ν ≠ ⌜seg⌟

xν), obeys an SIS dynamics (see S1 Text) with transmissibility λ, and can thus be computed as zwt = 1 − μ/λ. Secondly, when the whole set of segments circulates without wt (as compartment ⌜seg⌟), again the virus spreads as an SIS, this time with transmissibility (ρλ)v, and its endemic value can be predicted in the same fashion: zseg=1-ρ(ρλ)v. Interestingly, for high ρ, and a transmissibility λ > ρ−v/(v − 1), it turns out that zseg > zwt: the prevalence of the multipartite form is higher than that of the wt.

In order to fully characterize the endemic states, we resort now to stochastic spreading simulations (see Methods), focusing on the bipartite case (v = 2). A higher number of variants (v > 2) would not change the qualitatively picture; it would simply increase the possible values for the prevalence of hj by increasing the number of possible segments that survive through hijacking. We choose six points in the parameter space that lie in different phases (see Fig 2E), and for those values we carry out the simulations.

We firstly focus on homogeneous host population structures. The results are shown in Fig 3A. We characterize the endemic states in terms of their type (see Table 1), and plot their total prevalence, and the prevalence of the defective variants. In the points lying on the x-axis (labeled by wt) the defective variants have gone extinct, and the wt behaves like an SIS (states wt). The points lying on the diagonal have witnessed the extinction of the wt, and the defective variants are circulating together in the ⌜seg⌟ compartment (states seg). Their values match the theoretical prediction (dashed vertical lines). The solid vertical lines in Fig 3A are the theoretical predictions of wt prevalence. They match all equilibria of type both wt and hj, as in those cases the total prevalence coincides with the prevalence of the wt. The states all, whose total prevalence cannot be predicted analytically, have the highest prevalence. This picture further confirms the relationship between the theoretically predicted phases and the allowed endemic states (Fig 2D).

Previously we have stated that the critical surfaces are not sensitive to recovery rate and topology separately, but only to the parameter μ^ encoding both at the same time. Specifically, two populations with different recovery rate and contact heterogeneity, but with the same μ^, are indistinguishable from the point of view of their critical behavior. The endemic prevalences, however, break this symmetry, as one can see from Fig 3B, where we focus on P3 (Fig 2D) and get to μ^=0.25 both with one homogeneous (as in Fig 3B), and two heterogeneous population structures (with exponents γ = 3.5 and γ = 3.2). All the three configurations show all the equilibria, as expected by the critical behavior, but in the heterogeneous case the prevalence is consistently lower for each equilibrium.

Likelihood of different endemic states

Up to now, we have identified the phases (allowed endemic states) and computed the prevalence of such states. We now focus on the probability of occurrence of each state. For each of the usual points in Fig 2E, we show the probability of reaching each equilibrium in Fig 3C. This is achieved by counting the number of stochastic realizations that, starting from similar initial conditions, lead to that specific equilibrium (branching ratio of that equilibrium). Clearly, points P1 and P6 have only one endemic state, which then has a probability equal to one of being reached. For the other points, which have more than one possible endemic scenario, these probabilities are more informative, as they tell us the chances of the different viral forms taking over the population. We remark, however, that while both the critical behavior and the prevalence of the endemic states are inherent properties of the system that do not depend on the specific initial conditions chosen, this is not true for the probabilities of occurrence, which are clearly influenced by the initial infection status of the population. Given that, however, we computed them by seeding only one host in the ⌜all⌟ compartment to a susceptible population, we can say that our predictions are—at least qualitatively—reliable in an invasion scenario, in which the viral form is introduced by just one (or few) individuals.

Rise of multipartitism

Using the analytical characterization of the endemic phases and the numerical study of the equilibria, we now can investigate under which conditions the interplay between spreading dynamics and topology of contacts leads to the rise and persistence of multipartitism. We can also determine the nature of such emergence, in terms of a commensal relationship with the wt, or a true competitive advantage at the ecological level. For the sake of simplicity, we start by considering no differential transmissibility (ρ = 1): wt and segments have the same transmission probability. The relevant figures are Fig 2A and 2B, points P1, P2 and P3 in Figs 2E and 3. In this scenario, three phases are possible, and one crosses them all by increasing the transmissibility λ. The first one (λ just above T1) is the wt-phase, in which only the wt can circulate, and whenever a defective segment is produced, it quickly goes extinct. By increasing λ, we then encounter the contingent-phase. This phase predates the appearance of true multipartitism, as defective segments can hijack the wt to circulate. These segments cannot persist on their own, but the highly prevalent wt allows them to complete the replication cycle. At this stage, any defective segment is a commensal of wt, as the persistence of the former depends on the latter, while wt’s fitness remains unchanged. The emergence of multipartitism in this context is a contingent process: segments circulate simply because they are allowed to, causing no change to the overall fitness of the virus. Furthermore, there is no selective pressure towards complementation, as a combination of segments reconstructing the full genome without the presence of the wt (compartment ⌜seg⌟) would not be able to persist. This is confirmed by the functional form of T2 in Eq (4), which features no dependence on the number of complementing variants v: the survival of each mutant is independent of the presence of others, as effective replication and diffusion is wt-mediated. In other words, complementation would not make the variants fitter to the environment.

A further increase in λ takes us to the all-phase. Here, in addition to the commensal relationship between wt and segments, complementing variants are able to circulate on their own, without wt: the equilibrium seg emerges. Selection then imposes a bias on those segments that together reconstruct the genome, as they represent a new effective spreading configuration, and increase the overall viral fitness. They are thus advantaged with respect to purely commensal segments. This fitness-enhancing effect is quite straightforward: let us suppose that, due to viral or host population bottlenecks, or another stochastic event, wt prevalence goes down drastically, to the point where it is cleared from the system. In the contingent-phase this would lead to complete viral extinction, as all the segments would die out, too, as their persistence is linked to wt’s. In the all-phase, on the other hand, the virus is more resilient, as it can still circulate thanks to complementation. This time selective pressure toward complementation is well visible in the expression of Ts in Eq (4), which depends exponentially on the number of variants v. This fact is in qualitative agreement with results in, where it was shown that the larger the number of segments, the harder to reach the persistence of the segmented variants within a host. Even if the multipartite genome does not enjoy any microscopic advantage, it can rise to fixation if the monopartite virus undergoes stochastic extinction. Though fluctuations would also affect the multipartite form, and stochastic extinction of the monopartite form is not very likely, similar scenarios are relevant in virus evolution [45, 46] and cannot be discarded a priori.

The contingent-phase: A stepping stone towards multipartitism

The analysis of the prevalences (Fig 3A) confirms the evolutionary drivers behind the different phases, and adds information regarding crossed effects between viral types. Moreover, it allows us to uncover an evolutionary potential for multipartitism even in the absence of an explicit microscopic advantage. Let us focus on the contingent-phase (point P2 in Fig 3A), and the all state. The total prevalence of the latter state is higher than wt’s and hj’s in the same phase, implying that some hosts are infected by complementing segments without wt (compartment ⌜seg⌟), that is by a bona fide multipartite virus. Given that the multipartite form is not endemic in this phase—the contingent-phase relies on the wt for viral persistence—, this excess prevalence of the all state is a by-product of overall viral prevalence, rather than its driver: an extinction of the wt would quickly drive segmented variants to extinction. It is, however, an important one, as while it may not increase fitness in that specific environment, it permits the independent replication of the set of complementing variants; these are then able to invade other environments in which the wt could not persist, as we will see in the following.

Contact heterogeneity impairs multipartitism

Let us examine the effect of a heterogeneous contact network on the phases and equilibria above. As we have explained, in the phase space, topology is encoded in the parameter μ^. A low epidemic threshold is a well known feature of heterogeneous networks [41, 42, 44]. Specifically, power-law networks with γ < 3 exhibit a vanishing threshold as they grow larger, as the emergence of highly connected hubs ensures the persistence of the disease at any value of transmissibility, that is 〈k〉/〈k2〉 → 0 (and as a result, μ^→0) as the number of potential hosts grows, N → ∞. In our case this translates into T1, which is the epidemic threshold, going to zero for μ^→0. Also T2, Ts → 0. Further information is obtained when comparing their limit behaviors. As μ^ becomes smaller, T2/T1 increases but remains finite, while Ts/T1 → ∞. This implies that, the higher the heterogeneity of the network is, the more difficult it becomes for multipartitism to persist. Specifically, reaching the all-phase from the wt-phase would require an infinite relative increase (in the limit μ^→0) in transmissibility. Even when heterogeneity is not severe enough so as to cause the threshold to vanish, i.e., when γ > 3, heterogeneity makes it harder to sustain multipartitism, as both T2/T1 and Ts/T1 are decreasing functions of μ^.

Heterogeneity also modifies endemic prevalences and the branching ratios of equilibria. Let us examine point P3 (all-phase in Fig 2E): when the network is homogeneous, the highest branching ratio corresponds to all, and the equilibria containing segments together (hj) happens 20% of the time. When the network is heterogeneous, this fraction decreases, and wt quickly overtakes all in being the most probable outcome.

Summarizing, homogeneous contact patterns favor the emergence and persistence of multipartitism, while heterogeneous contacts hamper it. Qualitatively, this is the result of the complex interaction between the bottlenecks induced by between-host transmission and the presence of superspreaders, i.e., hosts that can potentially infect a large fraction of the population thanks to their high number of contacts. Combining this mechanism with a low MOI—and hence low λ—predicts an evolutionary radiation of multipartite viral forms linked to the rise and intensification of agricultural practices. In crops and cultivars, contacts among hosts are much more homogeneous (and often closer) than in wild settings, tremendously alleviating the requirements imposed by co-infection. Multipartite viruses adapted to the patchy distribution of wild hosts could have found it easy to propagate in regular, monospecific host populations which in all cases have closely related wild forms from which they departed through artificial selection.

Emergence of new multipartite phases through increased transmissibility of the segmented form

Up to this point, we have assumed that all viral forms have the same transmissibility and, still, fixation of the multipartite form cannot be fully discarded. Any additional advantage, however minor, of multipartitism will contribute to its ecological success, as we now discuss. We now set out to study the impact of enhanced transmissibility of defective variants, encoded in the parameter ρ being larger than one (ρ > 1). As ρ increases, the endemic state seg becomes more prevalent and more likely with respect to wt (see Fig 3C and 3D). Specifically, the value of the transmissibility for which zseg > zwt decreases as ρ increases, facilitating the predominance of the multipartite form (in Fig 3A point P3 has zwt < zseg, while P4 and P5 have zseg > zwt). Most importantly, ρ > 1 causes two new phases to emerge (Fig 3B, 3C and 3D), and both facilitate the rise of multipartitism by eliminating hj from their possible equilibria. One is the mix-phase, in which the virus circulates either as wt or as a multipartite. The mix-phase also presents an all endemic state that results from the interaction between the two former equilibria. Unlike in the all-phase, however, here the all equilibrium no longer indicates commensal relation. The second emerging phase is the seg-phase, in which only the complemented multipartite virus is able to circulate, while the monopartite version quickly goes extinct (yellow, and point P6 in Figs 2 and 3). This phase is of paramount importance because it lies in a parameter region where, without developing multipartitism, the virus would not be endemic. In addition, by prescribing the exclusive presence of the multipartite form, it allows to explain the phenomenology observed in nature, as the simultaneous presence of monopartite and defective-complementing forms of the same virus has been observed only in vitro [20, 48–50]. In vivo, viral species circulate as either pure monopartite (endemic state wt), or pure multipartite (seg). Furthermore, although in vitro several defective viral forms are generated and detected, and propagate along the wt (which would correspond to an in vitro

hj), this equilibrium is rarely found in wild plants. There is, however, an association between fully-fledged viruses and defective viral forms formally equivalent to the hj equilibrium: virus and viral satellites. Often, in addition, satellites modify the aetiology of viral infections, such that the transmissibility and the recovery rate might be affected by its presence in no particular direction, a phenomenon we do not consider in our model. There are other classes of hyperparasites that depend on a functional virus for replication (e.g. virophages or viroid-like satellites) whose ecological dynamics could, with appropriate modifications, be described in the framework discussed here. Interestingly, it has been proposed that virus-satellite associations, a typically unrelated tandem from a phylogenetic viewpoint, might evolve towards full co-dependence, and therefore be a possible, alternative evolutionary pathway to multipartitism.

Though our knowledge of existing viral forms is still incomplete and likely biased [16, 55], our results indicate that endemic states mixing monopartite and multipartite cognate forms (hj, all) need values of transmissibility difficult to sustain: endemicity could be achieved with lower values of transmissibility if the virus propagated only as a wild-type, while high values entail a cost that is usually compensated by decreasing infectivity. Albeit rare, however, these endemic states might act also as a stepping stone towards multipartitism even if they are only transiently present, as in the following example.

Consider a purely monopartite virus endemic in a plant population, in a specific environment, with parameters μ^A, λA as in point A in the inset of Fig 2C. μ^A is a combination of the recovery rate of the disease (characteristic of the host-virus interaction), and the between-plant contact network, driven by plant distribution and vector movements. A second population, occupying an adjacent geographic area, may have a different parameter (μ^B>μ^A), due to a different contact topology. As the inset of Fig 2C shows, the virus is able to colonize the second population only through an evolutionary process that increases its transmissibility up to at least λB′, so that point B′ is above the epidemic threshold (green path in the figure). It is reasonable to assume that the larger the increase in transmissibility required, the less likely this process is, given that the required mutation(s) are less likely and possibly more costly to maintain. The emergence of multipartitism decreases the evolutionary distance between the two states, increasing adaptability (magenta path in the figure). Random mutations, in fact, need to increase transmissibility from λA (wt-phase) to λB < λB′ (all-phase), where a complementing, multipartite version of the virus can emerge. Invasion of the second population is now possible, because the new viral forms effectively lowers the epidemic threshold in μ^B, thanks to the emergence of the seg-phase (point B). This simplistic example not only shows that multipartitism can emerge as a fitness-enhancing feature, but also that coexistence of monopartite and multipartite forms is a key stage in the evolutionary process, albeit possibly transient and short-lived. In addition to outlining the adaptive potential of multipartitism, this example elucidates the hampering effect of network heterogeneity. By increasing the distance between the wt-phase and the all-phase, network heterogeneity reduces the ratio λB′/λB, making multipartitism less advantageous. In conclusion, while making viral persistence overall easier, network heterogeneity curbs the potential of multipartitism as an effective adaptation strategy.

Conclusion

Multipartitism represents an example of a complex and as-of-today puzzling viral strategy. We have developed a framework that, starting from few key biological features, models the interaction between monopartite and multipartite forms, driven by the spreading dynamics on a host population. Despite assuming that multipartitism emerged from complementation between defective viral forms generated by the wt virus, as it has been observed in vitro, our results can be extended to other situations with relative ease. Most importantly, in addition, we have described how the structure of contacts among hosts drives the rise and persistence of the different viral forms. We have analytically characterized the parameter regions leading to viral persistence, in the form of wt only, of wt and defective segments, or segments only. We have also defined the different types of relationships between wt and segments, and specifically the presence or absence of selective pressure towards complementation, i.e., to witnessing the circulation of defective variants that may cooperatively reconstruct the whole genome. We have corroborated these findings through stochastic numerical simulations aimed at computing the prevalence of the different endemic states, and their probability of occurrence. As a result, we have been able to identify under which ecological conditions would multipartitism be a successful adaptive strategy, in the presence or absence of microscopic advantages, to new external conditions and environments characterized by variations in the topology of contacts between hosts. Defective particles generated through replication errors would start circulating by hijacking the wt. Subsequently, a complementing set of variants might form. Once that situation is achieved, even a small advantage in transmissibility (ρ > 1) would give an advantage to the multipartite form, which could anyway replace the monopartite form if chance causes the stochastic extinction of the latter. This sequence of events represents a plausible, parsimonious evolutionary pathway to the rise and persistence of multipartite viruses, and clarifies in which manner multipartitism might be an effective adaptive strategy at the ecological level.

We have also uncovered that while heterogeneous contact patterns among hosts favor viral persistence in general, they give a higher advantage to monopartite forms, by limiting the evolutionary and adaptative potential of multipartitism. Our model clearly lacks specific biological features that characterize different viruses, but that is a strength rather than a weakness, as it can be applied to a wide variety of settings with appropriate minimal modifications. We nonetheless explore additional realistic features in S1 Text, as nonhomogeneous recovery rates and nonindependent viral transmission.

Finally, it is worth discussing the effect of the interaction between the microscopic advantage and stochastic effects on multipartite fixation. The effect of the microscopic advantage, as quantified by our parameter ρ, becomes apparent in our current results, in terms of a much larger region of parameter space compatible with the fixation of the multipartite form (compare, for instance Fig 2B and 2C). Assuming a microscopic advantage therefore leads to a competitive advantage of multipartitism. A quantitative estimate of this competitive advantage, however, would require accounting for stochastic effects, an endeavor that goes beyond the current approach.

Despite not being able to formulate quantitative predictions, we are convinced that our framework provides an interesting qualitative picture of coexistence or substitution of different genomic architectures in a wide range of ecological environments. In this sense, we have uncovered evidence that the topology of contacts along which viruses spread may contribute to explaining why multipartite viruses preferentially infect plants. Our results lead us to conjecture that multipartite diversity and prevalence should have significantly increased together with the expansion of agriculture.

Materials and methods

Our goal is to derive the critical surfaces of Eqs (2)–(4) from the equation driving the dynamics of the system, namely Eq (1). We do that by starting from a simpler scenario, and incrementally adding features, up to the full model. Specifically, the first step consists in solving the model with no differential degradation (ρ = 1), no contact heterogeneity, and with only one segmented variant (v = 1). In the second step we generalize the result to a generic v, and in the third one we allow for heterogeneous contacts. In the last step we add differential degradation. A numerical validation of the critical surfaces is carried out in S1 Text.

⌜wt⌟ and ⌜all⌟ are the only infectious compartments, with prevalence x1 and x2, respectively. With neither differential degradation nor contact heterogeneity, Eq (1) reduces to

{x˙1=λ(1−x1−x2)x1+λ(1−λ)(1−x1−x2)x2−λx1x2−μx1x˙2=λ2(1−x1−x2)x2+λx1x2−μx2.

By summing these equations, we find that the equation for the total prevalence (z=defx1+x2) is z˙=λ(1-z)z-μz. This also follows from noticing that for v = 1 the total prevalence is also the total wt prevalence (see S1 Text). This means that the total prevalence behaves as a standard SIS, for which we know the epidemic threshold T1 = {λ = μ}, and the equilibrium above it. In addition, we know that just above T1 we are in the wt-phase. Hence, we have (zwt = 1 − μ/λ, x2,eq = 0). Studying the stability of this equilibrium gives us T2. Since the equation for z decouples from x1 and x2, it is convenient to study the system in (z, x2). Studying the sign of the eigenvalues of the Jacobian matrix reduces to studying ∂x˙2/∂x2<0 calculated in the equilibrium. This gives T2 = {λ = 2μ/(1 + μ)}. The details of the calculation are reported in the S1 Text.

We now generalize the previous result to an arbitrary v, while still assuming that all hosts have the same contact rate, that we can set to one with no loss of generality. Eq (1) simplifies to

x˙ν=∑βσΛνβσxβxσ+∑βΓνβxβ(1−∑σxσ)−μxν,(5)

whose Jacobian matrix is

Jνβ=∂x˙ν∂xβ=∑σ(Λν(βσ)-Γνσ)xσ+Γνβ(1-∑σxσ)-μδνβ,(6)

where Λν(βσ) = Λνβσ + Λνσβ. Firstly, we note that for v > 1 the total prevalence, now defined as z = ∑ν

xν, no longer behaves like an SIS, due to the presence of the compartment ⌜seg⌟. Indeed, one can show that, when summing over ν in Eq (1), the terms with Λ cancel out, as they pertain to interaction exclusively among infectious compartments, which by definition cannot change the total prevalence, and so all the contributions must cancel out. This is not the case however for the terms with Γ, so that the final equation is z˙=(1-z)∑β(Γx)β-μz, which does not decouple from xν. Interestingly, despite this breaking of the SIS symmetry, which was crucial to solve the v = 1 model, we can still prove that the values of T1, T2 found for v = 1 generalize to an arbitrary number of variants. We start from the first critical surface (T1). We compute the Jacobian matrix of Eq (6) in the dfs, i.e., xβ = 0∀β. We get J(dfs) = Γ − μ. We now argue that Γ, and therefore the Jacobian matrix, is upper triangular, thanks to the specific ordering of the compartments that we introduced. Γνβ is the rate at which a susceptible becomes a ⌜ν⌟, upon contact with a ⌜β⌟. For this to happen (Γνβ > 0), ⌜β⌟ must contain at least all the viral species ⌜ν⌟ contains. Hence, either β = ν (diagonal term), or β > ν. By the same reasoning, the diagonal terms are Γββ = λφ(β), where φ(β) is the number of viral species in ⌜β⌟, e.g. φ(⌜wt⌟) = 1, φ(⌜all⌟) = v+ 1. From these considerations, the spectrum of J(dfs) is {λφ(β)−μ;∀β}. Keeping in mind that λ < 1, we recover the first critical point: T1 = {λ = μ}.

Just above T1, ⌜wt⌟ is the only compartment with prevalence different from zero, hence it behaves like a standard SIS. Thanks to that we can compute Eq (6) in the wt-phase, and its spectrum. From that we find that the second critical surface is the same as for v = 1. The details of the calculation are in S1 Text.

We now build on the previous results, by adding heterogeneous contact rates. We work in the widely-used degree-block approximation [41–44, 57], assuming the contacts among agents are represented by an annealed network in which we assign each node a degree sampled from a power-law distribution with exponent γ: pγ(k) = Cγ

k−γ, where Cγ is the normalization factor. As customary, we assume γ > 2, so that the average degree is defined in arbitrary large populations. In the framework of annealed networks it makes sense to interpret k as a discrete number; one could also interpret it as a (continuous) coupling potential (either choice does not change the result found). We now directly compute the Jacobian of Eq (1), reported in Eq. (S.16) of S1 Text. The Jacobian is a matrix acting on a space which is the tensor product of the space of compartments, spanned by the Greek indices, and the space of degrees, spanned by the Latin indices. We can study its spectral properties on each space separately, using the previous results for the space of compartments. The full derivation is reported in the S1 Text.

Differential degradation ρ > 1 changes the matrices Γ, Λ, as reported in the S1 Text. The derivation is then similar to the case with ρ = 1.

Our model assumes that the transmission probability of one variant does not depend on the coinfecting variants. In reality, however, the number of viral particles a cell can produce in time is limited, and they are known to often spread in superinfection units. In S1 Text we investigate these aspects using a simple assumptions. We show that despite altering the specific values of the critical surfaces, they do not impact the qualitative behavior of the model.

Data in Fig 3 are produced through stochastic spreading simulations. Starting with a population of N = 6000, we infected the hosts with wt and both segments (⌜all⌟ for v = 2), and let the virus spread. We used an adaptation of the Gillespie algorithm [58, 59], to model both contacts among hosts, and contagion and recovery events. For each parameter configuration, we carried out 5000 simulations and kept only those reaching an endemic state other than the disease-free state, in order to discard instances of stochastic extinction, and focus only on the metastable equilibria which represent the attractors of the equations. We then used those simulations to compute prevalences and occurrence probabilities.