Red Queen dynamics in multi-host and multi-parasite interaction system

In host-parasite systems, dominant host types are expected to be eventually replaced by other hosts due to the elevated potency of their specific parasites. This leads to changes in the abundance of both hosts and parasites exhibiting cycles of alternating dominance called Red Queen dynamics. Host-parasite models with less than three hosts and parasites have been demonstrated to exhibit Red Queen cycles, but natural host-parasite interactions typically involve many host and parasite types resulting in an intractable system with many parameters. Here we present numerical simulations of Red Queen dynamics with more than ten hosts and specialist parasites under the condition of no super-host nor super-parasite. The parameter region where the Red Queen cycles arise contracts as the number of interacting host and parasite types increases. The interplay between inter-host competition and parasite infectivity influences the condition for the Red Queen dynamics. Relatively large host carrying capacity and intermediate rates of parasite mortality result in never-ending cycles of dominant types.

(Eq. S4) The population of a host type increases according to the effective growth rate (Gi) and decreases with the parasitic infections (functional response). The population of a parasite decreases by a constant death (dj) and increases with parasitic utilization of hosts (numerical response). We consider the following simplifying assumptions in the deterministic ODE system to investigate the basic dynamics of the Red Queen:  We assume equal number of host and parasite types, m=n, so that each host can have a unique specialist parasite.
 We assume equal host gross growth rates (ri=r for all i) and equal parasite death rates (dj=d for all j). The host types have the same characteristics (such as growth rate) but differ in susceptibility to parasites. Each parasite also shares the same characteristics (such as death rate) but targets different principal host. Note that the Red Queen dynamics still arise when these assumptions are relaxed.
 We assume logistic growth for host populations with carrying capacity K. Let We assume βii=1 and βik=0 for all k≠i, which implies that the densities of other hosts are marginal to alter the functional response curve for infecting host i. Red Queen is still possible when this assumption is relaxed (see Supplementary Fig. 5 where the value of βik is random with average αik).
 The parasitism efficiency matrix A=[aij] is diagonally dominant and symmetric, where Specifically, we set αij=α for i≠j and αij=α+1-αn for i=j. We hypothesize that differential susceptibility of hosts (infectivity of parasites) results in dominance replacement under this condition.
 The initial values are set to Hi(0)=H+0.001(i-1) for all i and Pj(0)=P+0.001(j-1) for all j to avoid unstable symmetric behavior. Without losing essential qualitative dynamics, we use H=P=0.01 in all the figures in the manuscript.
We define the Red Queen dynamics (or Red Queen cycles) as the case where all oscillating host/parasite densities exhibit dominance replacement by another host/parasite such that the amplitude of the oscillations are qualitatively identical but out-of-phase. Every host/parasite has the opportunity to be the sole dominant for a certain period of time. On the other hand, we generally refer to other types of oscillatory behavior as non-Red Queen dynamics.
Notes about the functional response: Type-I functional response denotes linear parasitic utilization of hosts. Type-II denotes parasitic utilization of hosts following a hyperbolic function. Type-III denotes parasitic utilization of hosts following a sigmoidal (S-curve) function. The curves of type-II and type-III converge to a saturation level; however, type-III includes parasite learning.

Numerical method
The system of ODEs is solved using Runge-Kutta 4 with stepsize=0.01. We use Berkeley Madonna (www.berkeleymadonna.com) to carry-out the computations. Maximum simulation time is t=30,000.

Details of the intensive simulation
In our search for the conditions that generate the Red Queen dynamics under the simplifying assumptions discussed above, we consider the following range of parameter values: where randN(μ,σμ) is a normal random number with mean μ and standard deviation σμ.

Additional captions for Figures 1 to 3
Parameter values for figures in the main text