Ecological model
We start with an ecological model of resident host-pathogen dynamics that assumes these populations are, respectively, genetically homogeneous. The ecological model underlies the evolutionary model we develop later. A complete description of the model, and the methods of analysis that follow, can be found in theSupplementary Information.
We consider a population of hosts classified according to their sex and disease status. At time t, there are Si=Si(t) sex-i individuals not infected by the pathogen, but susceptible to future infection (i=f for females, i=m for males). At time t there are also Ii=Ii(t) sex-i individuals who are not only infected with the pathogen but also able to transmit their infection to others. Our specific goal in this section is to develop a mathematical description of how the numbers of hosts in the various classes change over time.
The number of hosts in the population changes as a result of birth events. Following previous work44,45, we model the host mating rate using the harmonic mean of the population sizes of females and males. Assuming a one-to-one birth sex ratio, then newly born hosts of either sex join the population at rate (frac{b({S}_{f}+{I}_{f})({S}_{m}+{I}_{m})}{N}) where b>0, and N=N(t) denotes the total population size at time t. We assume that newborns produced by susceptible mothers are, themselves, susceptible. By contrast, we suppose that newborns produced by infected mothers acquire their mothers infection with probability v, where v is what we have called the vertical transmission rate31. Host number also changes because of death events. Hosts in every class experience natural mortality at per-capita rate N, where is a positive constant. Hosts infected by the pathogen also experience disease-related mortality at per-capita rate i (a measure of pathogen virulence) (Fig.7).
This model incorporates two sexes (females in red and males in blue) and vertical transmission (dashed line). The flow between compartments is represented by arrows and expressions next to each arrow represent the flow rate. Evolving phenotypes and drivers of their evolution are indicated in green and purple, respectively. Source data are provided as a Source Data file.
Numbers of hosts in any particular class changes as their disease-status changes. For example, we allow infected individuals to recover at per-capita rate i (a measure of host immunocompetence). We assume that, upon recovery, hosts move immediately into the appropriate susceptible group. In this way, we ignore the possibility that recovery implies immunity to subsequent infection. The disease status of hosts can also change because of horizontal disease-transmission events. We approach horizontal transmission in a standard way and assume that susceptible sex-i hosts acquire the pathogen horizontally from their infected sex-j counterparts at a total rate of SiijIj. Here, ij is a constant that reflects the transmissibility of the pathogen. We assume that when a host acquires an infection horizontally, it immediately becomes infectious (Fig.7).
The model described above is summarised mathematically using the following system of differential equations:
$$frac{d{S}_{f}}{dt}=frac{b({S}_{f}+(1-v){I}_{f})({S}_{m}+{I}_{m})}{N}+{gamma }_{f}{I}_{f}-{S}_{f}{beta }_{ff}{I}_{f}-{S}_{f}{beta }_{fm}{I}_{m}-mu N{S}_{f}$$
(1a)
$$frac{d{S}_{m}}{dt}=frac{b({S}_{f}+(1-v){I}_{f})({S}_{m}+{I}_{m})}{N}+{gamma }_{m}{I}_{m}-{S}_{m}{beta }_{mf}{I}_{f}-{S}_{m}{beta }_{mm}{I}_{m}-mu N{S}_{m}$$
(1b)
$$frac{d{I}_{f}}{dt}=frac{bv{I}_{f}({S}_{m}+{I}_{m})}{N}+{S}_{f}{beta }_{ff}{I}_{f}+{S}_{f}{beta }_{fm}{I}_{m}-({gamma }_{f}+{alpha }_{f}+mu N){I}_{f}$$
(1c)
$$frac{d{I}_{m}}{dt}=frac{bv{I}_{f}({S}_{m}+{I}_{m})}{N}+{S}_{m}{beta }_{mf}{I}_{f}+{S}_{m}{beta }_{mm}{I}_{m}-({gamma }_{m}+{alpha }_{m}+mu N){I}_{m}.$$
(1d)
Under a reasonable set of conditions, the previous system tends, over time, to an equilibrium state in which infections are endemic.
To study how pathogens disease-induced mortality and the hosts immune system respond to selection, we assume that each faces a life-history trade-off.
First, the pathogens ability to transmit horizontally trades off against the duration of any given infection it establishes. Following the previous authors30,46,47, we capture this trade-off by assuming
$${beta }_{ij}=beta ({alpha }_{j})=frac{{beta }_{max }{alpha }_{j}}{{alpha }_{j}+d}quad ,{{mbox{for}}},,j=f,;m,$$
(2)
where ({beta }_{max },,d , > , 0) are constants. Equation (2) implies that the nature of the trade-off faced by a pathogen is the same in both female and male hosts. Specifically, a pathogen can only increase its rate of horizontal transmission by increasing the disease-induced mortality rate experienced by its host (which, in turn, reduces the duration of infection). Equation (2) also says the horizontal transmission rate saturates at ({beta }_{max }) (independent of host sex), and does so more quickly as the parameter d is reduced (again, independent of host sex). Note also that Equation (2) does not depend on i: the sex of the susceptible host to whom the pathogen is transmitted.
For their part, hosts face a trade-off between investing resources in their immune system and their reproductive success. Increased immune investment is reflected in an increased recovery rate. To capture the hosts trade-off, then, we treat birth rate b as a decreasing function of the recovery rate. Moreover, we assume that the decrease in b is experienced by the host regardless of its disease status. In other words, we assume that cost associated with the immune system is an ongoing one, incurred mainly because of maintenance27 (this assumption model innate immunocompetence best) rather than being due to the activation that follows an infection48 (this assumption would model adaptive immunocompetence best). As noted in the Discussion, we relax this assumption in theSupplemental Material and compare the results for maintenance and activation costs. As an example, here, we point to evidence that shows female sex hormones enhance the immune system but simultaneously reduce the likelihood of conception and increase the chances of spontaneous abortion49,50,51. In mathematical terms, we capture the hosts trade-off using
$$b=b({gamma }_{f},{gamma }_{m})={b}_{max },{e}^{-{c}_{f}{gamma }_{f}^{2}},{e}^{-{c}_{m}{gamma }_{m}^{2}}$$
(3)
where ci reflects the rate at which fertility is reduced as sex-i immune function is increased (cost of recovery above). Equation (3) generalises the birth rate functions used previously27,48 to our sex-specific setting. The fact that b in this equation depends on both f and m reflects the fact that the reduced fertility of one mate affects the fertility of its partner16.
Our approach to modelling the co-evolution of host and pathogen is rooted in the adaptive-dynamics methodology52,53,54. For the pathogen population, we build a fitness expression that measures the success of a rare mutant strain in a population close to the endemic equilibrium established by the system (1) (indicated as ({bar{S}}_{i}), ({bar{I}}_{i}), and (bar{N})). Assuming that the mutant strain of pathogen is associated with a disease-induced mortality rate equal to ({tilde{alpha }}_{i}) in sex-i hosts, the number of mutant infections, ({tilde{I}}_{i}={tilde{I}}_{i}(t)) changes according to
$$frac{d{tilde{I}}_{f}}{dt}=frac{bv{tilde{I}}_{f}({bar{S}}_{m}+{bar{I}}_{m})}{bar{N}}+{bar{S}}_{f}beta ({tilde{alpha }}_{f}){tilde{I}}_{f}+{bar{S}}_{f}beta ({tilde{alpha }}_{m}){tilde{I}}_{m}-({gamma }_{f}+{tilde{alpha }}_{f}+mu bar{N}){tilde{I}}_{f}$$
(4a)
$$frac{d{tilde{I}}_{m}}{dt}=frac{bv{tilde{I}}_{f}({bar{S}}_{m}+{bar{I}}_{m})}{bar{N}}+{bar{S}}_{m}beta ({tilde{alpha }}_{f}){tilde{I}}_{f}+{bar{S}}_{m}beta ({tilde{alpha }}_{m}){tilde{I}}_{m}-({gamma }_{m}+{tilde{alpha }}_{m}+mu bar{N}){tilde{I}}_{m}.$$
(4b)
The system in (4) is linear and its long-term behaviour is determined by a dominant Lyapunov exponent of the mapping. We capture the information provided by the dominant Lyapunov exponent with the pathogen-fitness function, ({W}_{alpha }({tilde{alpha }}_{f},{tilde{alpha }}_{m},{alpha }_{f},{alpha }_{m})) using techniques laid out by the ref. 55 (see alsoSupplemental Information). When this function is greater than 1 the mutant invades and eventually displaces56 the resident strain associated with the i phenotype. When the function ({W}_{alpha }({tilde{alpha }}_{f},{tilde{alpha }}_{m},{alpha }_{f},{alpha }_{m})) is less than 1 the mutant does not invade and is eliminated from the population. With these facts in mind, we say that selection acts to move i in the direction given by the sign of (frac{partial {W}_{alpha }}{partial {tilde{alpha }}_{i}}{left|right.}_{tilde{alpha }=alpha }) where (tilde{alpha }=alpha) is shorthand for ({tilde{alpha }}_{i}={alpha }_{i}) for all i. Specifically, when this partial derivative is positive i is increasing, and when it is negative i is decreasing.
We follow a similar procedure for the host population by introducing, into the equilibrium population, a rare mutant-type host genotype that results in a recovery rate of ({hat{gamma }}_{i}) when expressed by sex-i hosts. We denote the numbers of susceptible and infected sex-i mutant-type hosts as ({hat{S}}_{i}) and ({hat{I}}_{i}), respectively. We assume that hosts are diploid, and so, strictly speaking, the hosts who contribute to ({hat{S}}_{i}) and ({hat{I}}_{i}) categories are heterozygotes (the numbers of homozygote mutants are negligible). While it remains rare, the dynamics of the mutant-host lineage can be described using
$$frac{d{hat{S}}_{f}}{dt}= frac{frac{b({hat{gamma }}_{f},{gamma }_{m})}{2}({hat{S}}_{f}+(1-v){hat{I}}_{f})({bar{S}}_{m}+{bar{I}}_{m})+frac{b({gamma }_{f},{hat{gamma }}_{m})}{2}({bar{S}}_{f}+(1-v){bar{I}}_{f})({hat{S}}_{m}+{hat{I}}_{m})}{bar{N}}\ +{hat{gamma }}_{f}{hat{I}}_{f}-{hat{S}}_{f}{beta }_{ff}{bar{I}}_{f}-{hat{S}}_{f}{beta }_{fm}{bar{I}}_{m}-mu bar{N}{hat{S}}_{f}$$
(5a)
$$frac{d{hat{I}}_{f}}{dt}= frac{frac{b({hat{gamma }}_{f},{gamma }_{m})}{2}v{hat{I}}_{f}({bar{S}}_{m}+{bar{I}}_{m})+frac{b({gamma }_{f},{hat{gamma }}_{m})}{2}v{bar{I}}_{f}({hat{S}}_{m}+{hat{I}}_{m})}{bar{N}}\ +{hat{S}}_{f}{beta }_{ff}{bar{I}}_{f}+{hat{S}}_{f}{beta }_{fm}{bar{I}}_{m}-({hat{gamma }}_{f}+{alpha }_{f}+mu bar{N}){hat{I}}_{f}$$
(5b)
$$frac{d{hat{S}}_{m}}{dt}= frac{frac{b({hat{gamma }}_{f},{gamma }_{m})}{2}({hat{S}}_{f}+(1-v){hat{I}}_{f})({bar{S}}_{m}+{bar{I}}_{m})+frac{b({gamma }_{f},{hat{gamma }}_{m})}{2}({bar{S}}_{f}+(1-v){bar{I}}_{f})({hat{S}}_{m}+{hat{I}}_{m})}{bar{N}}\ +{hat{gamma }}_{m}{hat{I}}_{m}-{hat{S}}_{m}{beta }_{mf}{bar{I}}_{f}-{hat{S}}_{m}{beta }_{mm}{bar{I}}_{m}-mu bar{N}{hat{S}}_{m}$$
(5c)
$$frac{d{hat{I}}_{m}}{dt}= frac{frac{b({hat{gamma }}_{f},{gamma }_{m})}{2}v{hat{I}}_{f}({bar{S}}_{m}+{bar{I}}_{m})+frac{b({gamma }_{f},{hat{gamma }}_{m})}{2}v{bar{I}}_{f}({hat{S}}_{m}+{hat{I}}_{m})}{bar{N}}\ +{hat{S}}_{m}{beta }_{mf}{bar{I}}_{f}+{hat{S}}_{m}{beta }_{mm}{bar{I}}_{m}-({hat{gamma }}_{m}+{alpha }_{m}+mu bar{N}){hat{I}}_{m}.$$
(5d)
The birth terms in the preceding system of equations reflect (a) the fact that the mutant host, while it is rare, mates only homozygous resident hosts and (b) only half of the matings between heterozygous mutants and homozygous residents result in mutant offspring. Since the dynamics described by (5) are linear, we can again measure fitness (this time for the host) using the dominant Lyapunov exponent. We summarise the relevant information contained in this exponent with the host fitness function ({W}_{gamma }({hat{gamma }}_{f},{hat{gamma }}_{m},{gamma }_{f},{gamma }_{m})), again using techniques outlined by ref. 55. In keeping with the description of pathogen evolution, we assert that the hosts i is increasing when (frac{partial {W}_{gamma }}{partial {gamma }_{i}}{left|right.}_{hat{gamma=gamma }}) is positive, and decreasing when this partial derivative is negative, where (hat{gamma }=gamma) is shorthand for ({hat{gamma }}_{i}={gamma }_{i}) for all i.
We want to identify where the action of selection takes the resident pathogen and host traits (i and i, respectively) in the long term. As mentioned above, the model is too complicated to support exact mathematical predictions. Consequently, our methods rely on numerical simulation implemented in Matlab57. All Matlab code is publicly available (see Code Availability).
The numerical simulation takes as its input a set of parameters and an initial estimate of the long-term result of selection on co-evolution of pathogen and host ({alpha }_{i}^{*}), and ({gamma }_{i}^{*}) for i=f, m. The estimate is updated by (i) finding the corresponding equilibrium solution to Equation (1) in a manner that verifies its asymptotic stability, (ii) using that equilibrium solution to estimate partial derivatives (frac{partial {W}_{alpha }}{partial {tilde{alpha }}_{i}}{left|right.}_{tilde{alpha=alpha }}) and (frac{partial {W}_{gamma }}{partial {hat{gamma }}_{i}}{left|right.}_{hat{gamma=gamma }}) for i=f, m, and finally (iii) incrementing or decrementing elements of the estimate following the sign of the appropriate partial derivative. Steps (i)(iii) are repeated until the absolute value of all partial derivatives is within a tolerance of zero. The result of the simulation is an estimate of the convergence stable58,59 co-evolutionary outcome, assuming f and m, and f and m can be adjusted independently. Importantly, this predicted co-evolutionary outcome also corresponds to a system in which the pathogen is established in a stable equilibrium population of hosts.
Finally, we verified numerically that the convergence-stable estimate corresponded to a two-dimensional evolutionarily stable result60 for pathogen and host, respectively. For this reason, we can also refer to predictions generated by our numerical simulation as a continuously stable state, in analogy to the definition established by ref. 61.
Further information on research design is available in theNature Research Reporting Summary linked to this article.
Continued here:
On maternity and the stronger immune response in women - Nature.com