abundance-density-tl.Rtex

%--------------------------------------------------------------------------------------------------------------------------------%
% Code and text for "A mechanistic model to compare the importance of interrelated population measures"
% by Tim CD Lucas, Hilde Herbots and Kate Jones
% See code chunk "settings" to see important switches (run large simulations or not etc.)
%
%---------------------------------------------------------------------------------------------------------------------------------%



\documentclass[11pt]{article}
\usepackage[sc]{mathpazo} %Like Palatino with extensive math support
\usepackage{fullpage}
\usepackage[authoryear,sectionbib,sort]{natbib}
\linespread{1.7}
\usepackage[utf8]{inputenc}
\usepackage{lineno}
\renewcommand{\refname}{Literature Cited}
\renewcommand{\cite}{\citep}
\usepackage[nomarkers, nolists,figuresonly]{endfloat}


\usepackage[width = 0.9\textwidth]{caption}
\usepackage{verbatim, graphicx, xcomment, array}
\usepackage{amsmath}

\DeclareGraphicsExtensions{pdf, eps, png}



\usepackage{booktabs}
\PassOptionsToPackage{hyphens}{url}
\usepackage[pdftex,hidelinks]{hyperref}



%%%%%%%%%%%%%%%%%%%%%
% Line numbering
%%%%%%%%%%%%%%%%%%%%%
%\usepackage{lineno}
% Please use line numbering with your initial submission and
% subsequent revisions. After acceptance, please turn line numbering
% off by adding percent signs to the lines %\usepackage{lineno} and
% to %\linenumbers{} and %\modulolinenumbers[3] below.

\title{Role of inter-related population-level host traits in determining pathogen richness and zoonotic risk}

% This version of the LaTeX template was last updated on
% April 28, 2017.

%%%%%%%%%%%%%%%%%%%%%
% Authorship
%%%%%%%%%%%%%%%%%%%%%
% Please remove authorship information while your paper is under review,
% unless you wish to waive your anonymity under double-blind review. You
% will need to add this information back in to your final files after
% acceptance.

\author{
  Tim C.D. Lucas\textsuperscript{1,2,$\ast$}\\
  Hilde M. Wilkinson-Herbots\textsuperscript{3}\\ 
  Kate E. Jones\textsuperscript{2,4,$\ast$}
}

\date{}



\begin{document}

\maketitle

\noindent{} 1. Oxford Big Data Institute, Li Ka Shing Centre for Health Information and Discovery, University of Oxford, Oxford, OX3~7FZ, United Kingdom.;

\noindent{} 2. Centre for Biodiversity and Environment Research, Department of Genetics, Evolution and Environment, University College London, Gower Street, London, WC1E~6BT, United Kingdom 

\noindent{} 3. Department of Statistical Science, University College London, Gower Street, London, WC1E~6BT, United Kingdom;

\noindent{} 4. Institute of Zoology, Zoological Society of London, Regent's Park, London, NW1~4RY, United Kingdom. 

\noindent{} \textbf{ORCIDs} Lucas, 0000-0003-4694-8107. Jones, 0000-0001-5231-3293.

\noindent{} $\ast$ Corresponding authors; e-mail: timcdlucas@gmail.com, kate.e.jones@ucl.ac.uk.

\bigskip

\textit{Word count}: 4371

\textit{Manuscript elements}: Figure~1, figure~2, table~1, online appendices~A and B (including figure~A1, figure~A2, figure~A3, figure~A4, table~A1, table~A2, table~A3 and table~A4). Figures~1 and  2 are to print in color.

\bigskip

\textit{Keywords}: Pathogen competition, zoonotic disease, metapopulations, host pathogen richness, bats, emerging infectious disease.

\bigskip

\textit{Manuscript type}: Article. %Or e-article, note, e-note, natural history miscellany, e-natural history miscellany, comment, reply, invited symposium, or countdown to 150.

\bigskip

\noindent{\footnotesize Prepared using the suggested \LaTeX{} template for \textit{Am.\ Nat.}}

\linenumbers{}
\modulolinenumbers[3]

\newpage{}

% running title: Population-level drivers of pathogen richness


%%begin.rcode settings, echo = FALSE, cache = FALSE, message = FALSE, results = 'hide'

####################################
### Important simulation options ###
####################################


# Compilation options
# Run simulations? This will take many hours
runAllSims <- FALSE

# Save raw simulation output
# This will take ~10GB or so.
# If false, summary statistics of each simulation are saved instead.
saveData <- FALSE

# Display extra figures in final documents?
extraFigs <- 'hide'


# How many cores do you want to use to run simulations?
nCores <- 8

##########################
### End options        ###
##########################


opts_chunk$set(cache.path = '.abundanceCache/', fig.path = 'figure/')

source('misc/KnitrOptions.R')


%%end.rcode


%%begin.rcode ggplot, cache = FALSE

source('misc/theme_tcdl.R')
theme_set(theme_grey() + theme_pub)

%%end.rcode

%%begin.rcode libs, cache = FALSE, result = FALSE

# My package. For running and analysing Epidemiological sims.
#   https://github.com/timcdlucas/metapopepi
library(MetapopEpi)

# Data manipulations
library(reshape2)
library(dplyr)

# Calc confidence intervals (could probably do with broom instead now.)
library(binom)

# To tidy up stats models/tests
library(broom)


# Run simulations in parallel
library(parallel)


# plotting
library(ggplot2)
library(palettetown)
# Get nice fonts in base graphics
library(showtext)

library(assertthat)

%%end.rcode

\section*{Abstract}

Zoonotic diseases are an increasingly important source of human infectious diseases, and host pathogen richness of reservoir host species is a critical driver of spill-over risk. 
Population-level traits of hosts such as population size, host density and geographic range size have all been shown to be important determinants of host pathogen richness. 
However, empirically identifying the independent influences of these traits has proven difficult as many of these traits directly depend on each other. 
Here we develop a mechanistic, metapopulation, susceptible-infected-recovered model to identify the independent influences of these population-level traits on the ability of a newly evolved pathogen to invade and persist in host populations in the presence of an endemic pathogen.
We use bats as a case study as they are highly social and an important source of zoonotic disease. 
We show that larger populations and group sizes had a greater influence on the chances of pathogen invasion and persistence than increased host density or the number of groups. 
As anthropogenic change affects these traits to different extents, this increased understanding of how traits independently determine pathogen richness will aid in predicting future zoonotic spill-over risk.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section*{Introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%\tmpsection*{General Intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A basic introduction to the field,
% comprehensible to a scientist in any discipline.

%1. Zoonotics important - determinants of pathogen richness important



% Pop-level is important as shown by corr studies.
% Sometimes the definitions are difficults
% While density been studied, population abundance has not despite being a fundementall characteristic.

Zoonotic diseases are a major source of human infectious disease \cite{jones2008global, taylor2001risk}.
Epidemics of emerging, zoonotic diseases pose a major threat to human health and economic development \cite{world2016ebola, ebola2015worldbank}.
The probability of zoonotic spill-over depends on, amongst other factors, the number of pathogen species carried by reservoir host species (pathogen richness) \cite{wolfe2000deforestation}. 
Empirical, comparative studies across reservoir host species, suggest that host morphological and life-history traits, such as large body size and longevity, correlate strongly with high pathogen richness \cite{kamiya2014determines, luis2013comparison}. 
However, traits related to reservoir host population biology are also expected to affect disease dynamics and therefore influence pathogen richness. 
Population-level traits such as increased host density \cite{morand1998density, kamiya2014determines, nunn2003comparative}, large geographic range size \cite{nunn2003comparative, turmelle2009correlates, kamiya2014determines} and greater population structure (nonrandom interactions between individuals) \cite{maganga2014bat, turmelle2009correlates} have been shown to correlate with high pathogen richness, although the evidence for a relationship with group size (number of individuals in a social group) has been equivocal in many studies \cite{vitone2004body, gay2014parasite, rifkin2012animals, nunn2003comparative, turmelle2009correlates}. 
Population size (total number of individuals), an important population-level trait, has rarely been included in comparative studies, despite its importance in describing epidemiological populations \cite{begon2002clarification}. 


%4. then talk about the limitations of taking a comparative approach 

Collinearity between explanatory variables is a common problem in correlative studies, and this issue is exacerbated when there are clear, causal relationships between explanatory variables.
There are two particularly clear relationships between the population-level traits associated with pathogen richness.
Firstly, host density, $d$, host population size, $N$, and geographic range size, $a$, are, by definition, linked by $d = N / a$ (see electronic supplementary material, table~S1 for all parameters used) and this relationship has broad empirical support \cite{blackburn2006variations}.
Secondly, host population size can be decomposed into two components, the number of groups, $m$, and the average size of a group, $n$, with $N = mn$.
%While less clear, geographic range size and host population structure are also related.
Correlative, comparative studies would be especially poor at identifying which, if any, of these traits causally affect pathogen richness.
This lack of discriminatory power is particularly important with respect to global change and its effects on zoonotic disease emergence.
Population-level traits such as population size and geographic range size, although interrelated, will respond differently to global change and the response will be species specific.
Some host species may suffer geographic range contractions while their density remains constant \cite{thomas2004extinction}.
Other host species might retain their geographic range but have a depressed population density \cite{craigie2010large}.
Only by knowing which of these interrelated traits control pathogen richness will we be able to predict future changes in pathogen richness.


%The alternative is mechanistic models.
Mechanistic models provide one method for comparing the importance of intrinsically related traits and can provide a deeper understanding of the system than correlative approaches.
Theoretical studies have established that a number of host population-level traits are important for epidemiological dynamics and the maintenance of pathogen richness.
In particular, host density, population structure and group size are well established as having central roles in pathogen dynamics \cite{colizza2007invasion, may1979population, anderson1979population, colizza2007invasion}.
A number of studies have found that increased host population structure can promote pathogen coexistence \cite{qiu2013vector, allen2004sis, nunes2006localized}.
While these studies have examined whether these population-level traits can promote pathogen richness, none have attempted to distinguish which might be the most important.
Mechanistic models that try to disentangle the interplay between population-level traits including host density, population size, geographic range size, group size and the number of groups are critically needed.


%7. what you are going to do that is new (ie tease apart these interdependencies) - using case study on bats (and why this is a good group to ask this particular question on)


Here, we use multipathogen models to individually vary host population-level traits.
We examine a simple deterministic model to establish whether a newly evolved pathogen can invade into an unstructured population in the presence of strong competition from an endemic pathogen strain. 
We then examine a stochastic, metapopulation model that was parameterised to broadly mimic wild bat populations.
We used bats as a case study as group (colony) size is very variable between bat species and bat colonies are often very stable \cite{kerth2011bats, mccracken1981social, jones2009pantheria, kerth2008causes}.
Furthermore, bats are particularly relevant in the context of zoonotic disease as they are thought to be reservoirs for a number of recent, important outbreaks \cite{calisher2006bats, li2005bats}.
We examined how the interrelated population-level traits affect the ability of a newly evolved pathogen to invade and persist in a population.
We used these simulations to examine two specific hypotheses.
First, we investigated whether host population size or density more strongly promotes the invasion of a new pathogen.
Secondly, we investigated whether the invasion of a new pathogen is more strongly promoted by group size or the number of groups.
We found that population size has a much stronger effect on the invasion of a new pathogen than host density.
We also found that increasing population size by increasing group size promotes pathogen invasion much more than increasing population size by increasing the number of groups.




%%begin.rcode SimSetup

  # These apply to both topo and disp sims. And probably should apply to extinction sims if I include them.
  # How long should the simulation last?
  nEvent <- 1.4e6
  # When should the invading pathogen be added.
  invadeT <- 6e5

  # Total time 
  TotalTime <- 75
  
  # Sims per parameter set.
  each <- 100
  
  # How often should the population be sampled. Only sampled populations are saved.
  sample <- 1000

  source('abundance-density-functions-tl.R')

%%end.rcode




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Constant density size. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%begin.rcode DensSimsFuncs

  #################################
  # Density sim definitions       #
  #################################

  # Keep density constant. Alter population size via colonySize

  # How many simulations to run?
  nDensSims <- 15 * each

  # Parameters for constant density
  #   Area and meanColonySize are the actal input arguments
  pop <- rep(c(2000, 4000, 8000, 16000, 32000), each = 3 * each)
  colonySize <- pop / 20
  
  # make a copy of colonySize for use in text.
  dep1 <- unique(colonySize)

  area <- pop / 0.8
  side <- sqrt(area)

  tran <- rep(c(0.1, 0.2, 0.3), times = 5 * each)

%%end.rcode

%%begin.rcode runDensSim, eval = runAllSims, cache = TRUE

# Create and set seed (seed object is used to set seed in each separate simulation.'
seed <- 1
set.seed(seed)

# If we want to save the data, make a directory for it.
if(saveData){
  dir.create('Data/fullsimoutput/')
}

# Run sims.
z <- mclapply(1:nDensSims, . %>% fullSim1, mc.preschedule = FALSE, mc.cores = nCores)

z <- do.call(rbind, z)

# Save summary data.
write.csv(z, file = 'Data/DensSims.csv')

%%end.rcode


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Constant density 2. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%begin.rcode Dens2SimsFuncs

  #################################
  # Density2 sim definitions      #
  #################################

  # Keep density constant. Alter population size via colonyNumber


  # How many simulations to run?
  nDens2Sims <- 15 * each

  # Parameters for constant density
  #   Area and meanColonySize are the actal input arguments
  pop <- rep(c(2000, 4000, 8000, 16000, 32000), each = 3 * each)
  colonySize <- 400
  colonyNumber <- pop / colonySize
  
  # make a copy of colony number for use in text
  dep2 <- unique(colonyNumber )

  area <- pop / 0.8
  side <- sqrt(area)




  tran <- rep(c(0.1, 0.2, 0.3), times = 5 * each)

%%end.rcode

%%begin.rcode runDens2Sim, eval = runAllSims, cache = TRUE

# Create and set seed (seed object is used to set seed in each separate simulation.'
seed <- 2
set.seed(seed)


# Run sims.
z <- mclapply(1:nDens2Sims, . %>% fullSim2, mc.preschedule = FALSE, mc.cores = nCores)

z <- do.call(rbind, z)

# Save summary data.
write.csv(z, file = 'Data/Dens2Sims.csv')

%%end.rcode




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Constant Population. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%begin.rcode PopSimsFuncs

  #################################
  # Population sim definitions    #
  #################################

  # Keep population constant. Alter density by area

  # How many simulations to run?
  nPopSims <- 15 * each

  # Parameters for constant density
  #   Area and meanColonySize are the actal input arguments
  pop <- 8000

  colonySize <- pop / 20

  area <- rep(c(40000, 20000, 10000, 5000, 2500), each = 3 * each)
  side <- sqrt(area)
  
  
  # save density values for use in text.
  dep3 <- c(40000, 20000, 10000, 5000, 2500)
  depDens <- pop / dep3

  tran <- rep(c(0.1, 0.2, 0.3), times = 5 * each)

%%end.rcode

%%begin.rcode runPopSim, eval = runAllSims, cache = TRUE

# Create and set seed (seed object is used to set seed in each separate simulation.'
seed <- 3
set.seed(seed)


# Run sims.
z <- mclapply(1:nPopSims, . %>% fullSim3, mc.preschedule = FALSE, mc.cores = nCores)

z <- do.call(rbind, z)

# Save summary data.
write.csv(z, file = 'Data/PopSims.csv')

%%end.rcode



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section*{Methods}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{figure}[t]
\centering
  \includegraphics[width=0.4\textwidth]{figure/SIRoption1}
  \caption[Schematic of the two-pathogen SIR model used]{
  Schematic of the two-pathogen SIR model used. 
  Individuals are in one of five epidemiological classes, susceptible (orange, $S$), infected with Pathogen~1, Pathogen~2 or both (blue, $I_1, I_2, I_{12}$, respectively) or recovered and immune from further infection (green, $R$).
  Transitions between classes occur as indicated by solid arrows and depend on transmission rate ($\beta$), coinfection adjustment factor ($\alpha$) and recovery rate ($\gamma$).
  Births ($\Lambda$)  and deaths ($\mu$) are indicated by dashed arrows.
  %Parameter symbols for transitions are indicated.
  Note that individuals in $I_{12}$ move into $R$, not back to $I_1$ or $I_2$. 
  That is, recovery from one pathogen causes immediate recovery from the other pathogen.
  }
\label{f:sir}
\end{figure}




\subsection*{Two pathogen SIR model}


We developed a multipathogen, susceptible-infected-recovered (SIR) compartment model to examine the probability that a newly evolved pathogen would invade and persist into a population in the presence of an identical, endemic pathogen.
Susceptible individuals were counted in class $S$ (figure~\ref{f:sir}) while infected individuals were counted in classes $I_1$, $I_2$ and $I_{12}$, being individuals infected with Pathogen~1, Pathogen~2 or both, respectively.
Recovered individuals, $R$, were immune to both pathogens, even if they had only been infected by one (i.e.\ there was complete cross-immunity).
Furthermore, recovery from one pathogen moved an individual into the recovered class, even if the individual was coinfected (figure~\ref{f:sir}).
Though our study was restricted to two pathogens, this modelling choice allows the model to be easily expanded to include many pathogens. 
Our assumption of immediate recovery from all other pathogens is likely to be reasonable \cite{munywoki2015influence} as any up-regulation of innate immune response or acquired immunity would affect both pathogens equally.
The coinfection rate (the rate at which an infected individual is infected with a second pathogen) was adjusted compared to the infection rate by a factor $\alpha$.
Birth rate ($\Lambda$) was set equal to the death rate ($\mu$), meaning the population did not systematically increase or decrease.
New born individuals entered the susceptible class.
Infection and coinfection were assumed not to cause increased mortality as bats show no clinical signs of infection for a number of viruses \cite{halpin2011pteropid, deThoisy2016bioecological}.


Population structure is present in bat populations as demonstrated by the existence of subspecies, measurements of genetic dissimilarity and behavioural studies  \cite{kerth2011bats, mccracken1981social, burns2014correlates}.
Therefore assuming a fully-mixed population is an oversimplification.
Consequently, the population was modelled as a metapopulation network with groups being nodes and dispersal between groups being indicated by edges (electronic supplementary material, figure~S1).
Individuals within a group interacted randomly so that the group was fully mixed.
Dispersal occurred at a rate $\xi$ between groups connected by an edge in the network.
The dispersal rate from a group $y$ with degree $k_y$ to group $x$ was $\xi / k_y$.
Note that this rate was not affected by the degree and size of group $x$ and the total rate of dispersal was not affected by the degree of a group.
We examined the full model using stochastic, continuous-time simulations, in \emph{R} \cite{R, metapopepi}.
The full details of the model are given in electronic supplementary material~S1.


\subsection*{Deterministic model}
\label{s:determ}

We examined a single-group, deterministic model to establish a baseline expectation for whether a newly evolved pathogen strain could invade into a population in the presence of an identical strain given the assumptions of our two-pathogen SIR model (for details see electronic supplementary material~S2).
If we first consider the endemic pathogen (Pathogen~1) we have a typical SIR model with vital dynamics \cite{Kermack1932} with equilibrium values $S^\ast = \frac{\mu + \gamma}{\beta}$ and $I_1^\ast = \frac{\Lambda n}{\gamma + \mu} - \frac{\mu }{\beta}$ where $\beta$ and $\gamma$ are the transmission and recovery rates and $n$ is the group size.
When Pathogen~2 is introduced, its rate of change can be written as

\begin{align}
\frac{dI_2}{dt} & = \beta S^\ast I_2 + \alpha \beta I_1^\ast I_2 - (\gamma + \mu) I_2 \label{path2I}
\end{align}
which is greater than zero when $\alpha\left( \Lambda R_0 - \mu \right) I_2 > 0$ (with $R_0 = \frac{\beta N }{\gamma + \mu}$ being the basic reproduction number and being equal for the two identical pathogens).
As $\Lambda = \mu$ due to the assumption of a stable population size,  $\Lambda R_0 - \mu$ is greater than zero,
Therefore, $\frac{dI_2}{dt}$ is greater than zero provided $\alpha$ is greater than zero.
That is, provided cross-immunity is not complete, Pathogen~2 will invade in this deterministic model.
This means that it is only stochastic extinction that would prevent a pathogen from invading into a population.
Our more detailed simulations are therefore examining how effectively different population-level traits alleviate stochastic extinction or allow longer term persistence.

\subsection*{Parameter selection}
\label{s:paramSelect}

While some fixed parameters were chosen to approximate those found in wild bat populations, others were chosen for modelling reasons. 
Assuming equal birth and death rates, $\Lambda$ and $\mu$, were both set to 0.05 per year giving a generation time of 20 years \cite{jones2009pantheria, kerth2008causes}.
Setting birth and death rates equal gives a stochastically varying population size but given the length of the simulations groups were very unlikely to go extinct.
Although it is difficult to directly estimate infection durations in wild bat populations \cite{plowright2015ecological}, evidence suggests these can be on the scale of days \cite{amengual2007temporal} up to months or years \cite{peel2012henipavirus, sohayati2011evidence, rahman2010characterization}. 
Here we set the infection duration, $\frac{1}{\gamma}$, to one year which is a long lasting, but non-chronic, infection. 
As estimating transmission rates is particularly difficult we used three values of the transmission rate, $\beta$:  0.1, 0.2 and 0.3.
These values were chosen as very high values of $R_0$ were required so that a reasonable number of simulations experienced an invasion of Pathogen~2.
The coinfection adjustment parameter, $\alpha$, was set to 0.1.
%Therefore, an individual with a single infection was 90\% less likely to gain a second infection.
The deterministic model showed that $\alpha = 0$ and $\alpha > 0$ are qualitatively different with the number of individuals infected with  Pathogen~2 being stable and increasing respectively.
The case where Pathogen~2 does not invade and spread ($\alpha = 0$) is not important for pathogen richness so we chose a small, non-zero value for $\alpha$.
The dispersal rate, $\xi$, was set to 0.01 which yields 17\% of individuals dispersing in their lifetime.
This relates to a species with juvenile dispersal of a proportion of males and very few females \cite{kerth2008causes, chen2008sex}.

%\subsection*{Population-level factors} % This isn't a full subsection...

The effect of geographic range size on disease dynamics occurred through changes in the metapopulation network.
Dispersal was only allowed to occur between two groups if they were connected nodes in the metapopulation network.
The metapopulation network was created for each simulation by randomly placing groups in a square space which represented the geographic range of the species (electronic supplementary material, figure~S1).
Groups within a threshold distance of each other were connected in the metapopulation network.
This meant that the metapopulation network was not necessarily connected (made up of a single connected component).
To ensure connected metapopulation networks would have required repeatedly resampling the group locations until a connected metapopulation  occurred.
However, this would have biased the mean degree, $\bar{k}$.
Therefore, it was considered preferential to keep the unconnected networks.



\subsection*{Experimental setup}

A total of 4500 simulations were run.
In each simulation, each group in the host population was seeded with 20 individuals infected with Pathogen~1. 
A `burn-in' of \rinline{invadeT} events was run to allow Pathogen~1 to spread and reach equilibrium. 
Plotting of preliminary simulations was used to determine that \rinline{invadeT} events was enough to ensure equilibrium.
After this burn-in, five host individuals infected with Pathogen~2 were added to one randomly selected group. 
The invasion and persistence of Pathogen~2 was considered successful if any individuals infected with Pathogen~2 remained at the end of the simulation.
As simulations with many individuals and infection events had more events per unit time, the end of the simulation had to be defined in terms of time rather than the number of events.
Simulations were run until \rinline{TotalTime} years had elapsed since the introduction of Pathogen~2.
This value was chosen so that pathogens had to persist for multiple host generations in order to be considered persistent.
%The results were also checked to confirm that decreasing this value did not qualitatively change the results. 

Three sets of 1500 simulations were run in which three population parameters were varied: (\emph{i})  geographic range size (with corresponding change in host density), (\emph{ii}) group size (with corresponding change in population size), and (\emph{iii})  the number of groups (with corresponding change in population size).
To allow comparisons between simulation sets, the parameter that was being varied in each set was assigned its default value multiplied by 0.25, 0.5, 1, 2 and 4. 
To examine whether host density  had a stronger effect on pathogen invasion than population size, results from simulation set \emph{i} were compared to those from sets \emph{ii} and \emph{iii}.
To examine whether group size or the number of groups was the more important component of population size, results from simulation set \emph{ii} were compared to those from \emph{iii}.


%The default geographic range size was 10000 km\textsuperscript{2} with a range of 2500--40000 km\textsuperscript{2}.
The spatial scale is arbitrary; only the ratio between the geographic range size and the threshold distance for groups being connected in the metapopulation network had any effect on simulation outcomes.
We parameterised the spatial scale by arbitrarily selecting a threshold distance of 100 km. .
The default (10000 km\textsuperscript{2}) and upper and lower bounds of the geographic range size (2500--40000 km\textsuperscript{2}) were then selected to maximise the range of $\bar{k}$ (electronic supplementary material, figure~S2) while not having many simulations with networks that were unconnected.
That is, we aimed for low host density populations to have relatively sparse population networks, while high host density populations had fully-connected metapopulation networks.
This reflects the existence of both isolation-by-distance \cite{chen2008sex, o2015genetic, vonhof2015range} and panmixia \cite{peel2013continent, petit1999male, calisher2006bats} in bat species.
The default group size was 400 with a range of 100--1600 which is representative of many bat species \cite{jones2009pantheria}.
The default number of groups was 20 with a range of 5--80.
This minimum value is close to the minimum possible for the population to still be considered a metapopulation, while the maximum value was constrained by computational costs.

In the first set of simulations (\emph{i}), host density was varied by keeping population size constant ($N$ = 8000, $n$ = 400, $m$ = 20) while varying geographic range size.
The values of geographic range size used were 2500, 5000, 10000, 20000 and 40000 km\textsuperscript{2} which gave density values of \rinline{depDens[5:2]} and \rinline{depDens[1]} animals per km\textsuperscript{2}.
In the second set of simulations (\emph{ii}), population size was varied by changing group size while the number of groups was kept constant.
To keep host density constant, geographic range size was increased as population size increased.
The values of group size used were 100, 200, 400, 800 and 1600 while geographic range size was set to 2500, 5000, 10000, 20000 and 40000 km\textsuperscript{2}.
This gave population sizes of 2000, 4000, 8000, 16000 and 32000 while host density remained at 0.8 hosts per km\textsuperscript{2}.
In the third set of simulations (\emph{iii}), population size was varied by changing the number of groups while group size was kept constant.
Again, to keep host density constant, geographic range size was increased as population size increased.
The numbers of groups used were 5, 10, 20, 40 and 80 while geographic range size was set to 2500, 5000, 10000, 20000 and 40000 km\textsuperscript{2}.
Again, this gave population sizes of 2000, 4000, 8000, 16000 and 32000 while host density remained at 0.8 hosts per km\textsuperscript{2}.








%%begin.rcode calcMeans
# Constant density, altered colony size
# Need to remove rows with errors to auto find classes of columns
dens1text <- readLines('Data/DensSims.csv')
rmerrors <- dens1text[!grepl('Error in', dens1text)]
rmerrors <- rmerrors[nchar(rmerrors) > 3]
dens1 <- read.csv(text = rmerrors) %>% 
           dplyr::select(-X) 

dens1Means <- dens1 %>%
                #dplyr::filter(nExtantDis > 0) %>%
                group_by(transmission, colonySize, colonyNumber) %>%
                summarise(success = sum(nExtantDis == 2), sampleSize = n(), pop = mean(pop)) %>%
                ungroup %>% 
                cbind(., vary = 'colonySize', binom.confint(.$success, .$sampleSize, conf.level = 0.95, methods = "exact"))



# fit glms to each transmission value

colonySizeModel   <- dens1 %>% 
                       mutate(invasion = nExtantDis == 2 ) %>%
                       mutate(ValueChange = log2(colonySize / median(colonySize))) %>%
                       #filter(nExtantDis != 0) %>%
                       group_by(transmission) %>%
                       do(glm(invasion ~ ValueChange, data = ., family = 'binomial') %>% tidy(conf.int = TRUE)) 

colonySizePredictions <- dens1 %>% 
                           mutate(invasion = nExtantDis == 2 ) %>%
                           mutate(ValueChange = colonySize / median(colonySize)) %>%
                           #filter(nExtantDis != 0) %>%
                           group_by(transmission) %>%
                           do(augment(glm(invasion ~ log2(ValueChange), data = ., family = 'binomial'),
                                      newdata = data.frame(ValueChange = seq(min(.$ValueChange), max(.$ValueChange), 
                                                                                 length.out = 1000)),
                                      type.predict = 'response')) %>%
                           mutate(vary = 'colonySize')




# Constant density, altered colony number

dens2 <- read.csv('Data/Dens2Sims.csv') %>% dplyr::select(-X) %>% as.data.frame

dens2Means <- dens2 %>%
                #dplyr::filter(nExtantDis > 0) %>%
                group_by(transmission, colonySize, colonyNumber) %>%
                summarise(success = sum(nExtantDis == 2), sampleSize = n(), pop = mean(pop)) %>% 
                ungroup %>%
                cbind(., vary = 'colonyNumber',  binom.confint(.$success, .$sampleSize, conf.level = 0.95, methods = "exact"))



# Fit glms to each transmission value

colonyNumberModel <- dens2 %>% 
                       mutate(invasion = nExtantDis == 2 ) %>%
                       mutate(ValueChange = log2(colonyNumber / median(colonyNumber))) %>%
                       #filter(nExtantDis != 0) %>%
                       group_by(transmission) %>%
                       do(glm(invasion ~ ValueChange, data = ., family = 'binomial') %>% tidy(conf.int = TRUE))
                       
                       
colonyNumberPredictions <- dens2 %>% 
                             mutate(invasion = nExtantDis == 2 ) %>%
                             mutate(ValueChange = colonyNumber / median(colonyNumber)) %>%
                             #filter(nExtantDis != 0) %>%
                             group_by(transmission) %>%
                             do(augment(glm(invasion ~ log2(ValueChange), data = ., family = 'binomial'),
                                        newdata = data.frame(ValueChange = seq(min(.$ValueChange), max(.$ValueChange), 
                                                                                   length.out = 1000)),
                                        type.predict = 'response')) %>%
                             mutate(vary = 'colonyNumber')

                       
# Constant population, altered area

pop <- read.csv('Data/PopSims.csv') %>% dplyr::select(-X) %>% as.data.frame

popMeans <- pop %>%
                #dplyr::filter(nExtantDis > 0) %>%
                group_by(transmission, dens) %>%
                summarise(success = sum(nExtantDis == 2), sampleSize = n()) %>%
                ungroup %>%
                cbind(., binom.confint(.$success, .$sampleSize, conf.level = 0.95, methods = "exact"))


# fit glms to each transmission value

areaModel <- pop %>% 
               mutate(invasion = nExtantDis == 2 ) %>%
               mutate(ValueChange = log2(dens / median(dens))) %>%
               #filter(nExtantDis != 0) %>%
               group_by(transmission) %>%
               do(glm(invasion ~ ValueChange, data = ., family = 'binomial') %>% tidy(conf.int = TRUE))


areaPredictions <- pop %>% 
                     mutate(invasion = nExtantDis == 2 ) %>%
                     mutate(ValueChange = dens / median(dens)) %>%
                     #filter(nExtantDis != 0) %>%
                     group_by(transmission) %>%
                     do(augment(glm(invasion ~ log2(ValueChange), data = ., family = 'binomial'),
                                newdata = data.frame(ValueChange = seq(min(.$ValueChange), max(.$ValueChange), 
                                                                           length.out = 1000)),
                                type.predict = 'response')) %>%
                     mutate(vary = 'Density')



%%end.rcode



%%begin.rcode descriptive

meanDefault <- rbind(dens1, dens2, pop) %>%
                 filter(colonyNumber == 20, colonySize == 400, dens == 0.8) %>%
                 group_by(transmission) %>%
                 summarise(mean = mean(nExtantDis == 2))






transTest <- rbind(dens1, dens2, pop) %>%
               filter(colonyNumber == 20, colonySize == 400, dens == 0.8) %>%
               mutate(invasion = nExtantDis == 2) %>%
               glm(invasion ~ transmission, data = .) %>%
               tidy



nAllExtinct <- rbind(dens1, dens2, pop) %>%
                 filter(nExtantDis == 0) %>%
                 nrow


whichTransAllExtinct <- rbind(dens1, dens2, pop) %>%
                          filter(nExtantDis == 0) %>%
                          group_by(transmission) %>%
                          summarise(n = n())

transExtinctTest <- rbind(dens1, dens2, pop) %>%
                      mutate(doubleExtinct = nExtantDis == 0) %>%
                      glm(doubleExtinct ~ transmission, data = .) %>%
                      tidy



smallAllExtinct <- rbind(dens1, dens2, pop) %>%
                          filter(nExtantDis == 0, colonySize == 100) %>%
                          summarise(n = n())

fewAllExtinct <- rbind(dens1, dens2, pop) %>%
                          filter(nExtantDis == 0, colonyNumber == 5) %>%
                          summarise(n = n())


%%end.rcode





%%begin.rcode plotKcapt

plotKcapt <- '
Change in average metapopulation network degree ($\\bar{k}$) with increasing geographic range size. 
Bars show the median, boxes show the interquartile range, vertical lines show the range and grey dots indicate outlier values.
Notches indicate the 95\\% confidence interval of the median.
All simulations had 20 groups, meaning 19 is the maximum value of $\\bar{k}$.
'

plotKtitle <- 'Change in average network degree with increasing geographic range size'
%%end.rcode

%%begin.rcode plotK, fig.cap = plotKcapt, fig.scap = plotKtitle, out.width = '0.4\\textwidth', cache = FALSE, plot.height = 2.3, fig.show = 'hide'

ggplot(dens1, aes(x = factor(area), y = meanK, colour = 'a', fill = 'a')) +
  geom_boxplot(outlier.size = 1, size = 0.3, outlier.colour = grey(0.3), 
    notch = TRUE, width = 0.4, notchwidth = 0.6) +
  scale_colour_manual(values = pokepal('vileplume')[4]) +
  scale_fill_manual(values = pokepal('vileplume')[7]) +
  stat_summary(geom = "crossbar", width = 0.2, fatten = 2, 
    fun.data = function(x){ return(c(y = median(x), ymin = median(x), ymax = median(x))) }) +
  theme(legend.position = 'none', panel.grid.major.x = element_blank()) + 
  scale_x_discrete(labels = c('2500', '5000', '10000', '20000', '40000')) +
  xlab(expression(paste('Area ', (km^2)))) +
  ylab(expression(paste('Mean degree, ', italic(bar(k))))) +
  ylim(0, 20)

%%end.rcode








%%begin.rcode transMeansCapt

transMeansCapt <- '
Comparison of the effect of host population size on probability of invasion when population size is altered by changing group size or number of groups.
Relationships are shown separately for each transmission value, $\\beta$.
It can be seen that changes in group size give a much greater increase in invasion probability than changes in group number.
Note that this is the same data as figure~\\ref{fig:plotValueChangeMeans} but with the $x$-axis scaled by population size, rather than relative parameter change.
'

transMeansTitle <- 'Comparison of the probability of the affect of group size and number of groups on probability of invasion'

%%end.rcode


%%begin.rcode plotTransMeans, fig.cap = transMeansCapt, fig.scap = transMeansTitle, out.width = '0.4\\textwidth', fig.show = TRUE, cache = FALSE, fig.height = 4, fig.show = 'hide'

# Convert to percentage change in density, colonysize or colony number.

d2 <- rbind(dens1Means %>% 
             mutate(ValueChange = colonySize/median(colonySize), vary = 'colonySize') %>%
             dplyr::select(pop, transmission, mean, lower, upper, vary),
           dens2Means %>% 
             mutate(ValueChange = colonyNumber/median(colonyNumber), vary = 'colonyNumber') %>%
             dplyr::select(pop, transmission, mean, lower, upper, vary)
     )

d2$vary <- factor(d2$vary, levels = unique(d2$vary)[c(2, 3, 1)])
d2$transmission <- factor(d2$transmission, labels = c('beta == 0.1', '0.2', '0.3'))

# Join predictions


colonySizePredictions2 <- dens1 %>% 
                             mutate(invasion = nExtantDis == 2 ) %>%
                             #filter(nExtantDis != 0) %>%
                             group_by(transmission) %>%
                             do(augment(glm(invasion ~ pop, data = ., family = 'binomial'),
                                        newdata = data.frame(pop = seq(min(.$pop), max(.$pop), 
                                                                                   length.out = 1000)),
                                        type.predict = 'response')) %>%
                             mutate(vary = 'colonySize')
                       
colonyNumberPredictions2 <- dens2 %>% 
                             mutate(invasion = nExtantDis == 2 ) %>%
                             #filter(nExtantDis != 0) %>%
                             group_by(transmission) %>%
                             do(augment(glm(invasion ~ pop, data = ., family = 'binomial'),
                                        newdata = data.frame(pop = seq(min(.$pop), max(.$pop), 
                                                                                   length.out = 1000)),
                                        type.predict = 'response')) %>%
                             mutate(vary = 'colonyNumber')

percChangePreds2 <- rbind(colonySizePredictions2, colonyNumberPredictions2)

percChangePreds2$transmission <- factor(percChangePreds2$transmission, labels = c('beta == 0.1', '0.2', '0.3'))


ggplot(d2, aes(x = pop, y = mean, colour = vary)) +
  geom_line(data = percChangePreds2, aes(x = pop, y = .fitted, colour = vary)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper),
              width=.2) +
  ylab('Prob. Invasion') +
  xlab('Population size') +
  ylim(0, 1) +
  scale_colour_manual(name = "Population-level trait", 
                      labels = c('Number of groups', 'Group size'), 
                      values = pokepal('Swampert', spread = 3)[2:1]) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c('0.0', '0.5', '1.0')) +
  scale_x_continuous(limits = c(0, 32000), 
                     breaks = c(0, 1e4, 2e4, 3e4), 
                     label = c('0', '10000', '20000', '30000')) +
  facet_grid(transmission ~ ., labeller = label_parsed) +
	theme(legend.position = 'bottom', legend.direction = 'vertical')

%%end.rcode








\subsection*{Statistical analysis}



For each set of simulations, we fitted binomial GLMs in R \cite{R} with the proportion of invasions of Pathogen~2 as the response variable.
To enable comparison between GLMs we divided the explanatory variables by their default values and $\log_2$ transformed.
The explanatory variables for all three sets of simulations therefore became evenly spaced between -2 and 2.
To investigate whether an increase in host population size created a stronger increase in probability of pathogen invasion than an equal increase in host density we compared the size (and 95\% confidence intervals) of the regression coefficients of host density (\emph{i}) to group size (\emph{ii}) and number of groups (\emph{iii}). 
To examine whether an increase in group size creates a stronger increase in invasion probability than a proportionally equal increase in number of groups we compared regression coefficients of group size (\emph{ii}) to number of groups (\emph{iii}).
%In \rinline{nAllExtinct} out of 4500 simulations, both pathogens went extinct  (electronic supplementary material, tables~S2--4).
%Of these simulations, 29 were in set \emph{ii} and eight were in set \emph{iii}.
%All \rinline{nAllExtinct} simulations had either the smallest group size (\rinline{smallAllExtinct} simulations in set \emph{ii} with $n$ = 100,) or the fewest number of groups (eight simulations in set \emph{iii} with $m = 5$).
%However, the number of simulations in which both pathogens went extinct did not significantly depend on transmission rate (GLM: coefficient = \rinline{transExtinctTest$estimate[2]}, $p$ = \rinline{transExtinctTest$p.value[2]}).
%Results from these simulations were retained in the statistical analysis as they still represented the case where Pathogen~2 did not invade and persist.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section*{Results}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%begin.rcode valueChangeMeansCapt

valueChangeMeansCapt <- '
Comparison of the effect of population-level traits on probability of invasion. 
Population-level traits are group size (green lines, squares), number of groups (blue lines, circles) and host density (yellow lines, triangles).
The $x$-axis shows the change ($\\times$0.25, 0.5, 1, 2 and 4) in each of these traits relative to the default value.
Default values are: number of groups = 20, group size = 400 and host density = 0.8 animals per km\\textsuperscript{2}.
Each point is the mean of 100 simulations and bars are 95\\% confidence intervals.
Each curve was obtained by fitting a binomial GLM.
Relationships are shown separately for each transmission value, $\\beta$.
'

valueChangeMeansTitle <- 'Comparison of the effect of group size, number of groups and host density on probability of invasion'

%%end.rcode







%%begin.rcode plotValueChangeMeans, fig.cap = valueChangeMeansCapt, fig.scap = valueChangeMeansTitle, out.width = '0.5\\textwidth', cache = FALSE, fig.height =  6

# Convert to percentage change in density, colonysize or colony number.

d <- rbind(popMeans %>% 
             mutate(ValueChange = dens/median(dens), vary = 'Density') %>%
             dplyr::select(ValueChange, transmission, mean, lower, upper, vary),
           dens1Means %>% 
             mutate(ValueChange = colonySize/median(colonySize), vary = 'colonySize') %>%
             dplyr::select(ValueChange, transmission, mean, lower, upper, vary),
           dens2Means %>% 
             mutate(ValueChange = colonyNumber/median(colonyNumber), vary = 'colonyNumber') %>%
             dplyr::select(ValueChange, transmission, mean, lower, upper, vary)
     )

d$vary <- factor(d$vary, levels = unique(d$vary)[c(2, 3, 1)])
d$transmission <- factor(d$transmission, labels = c('beta == 0.1', '0.2', '0.3'))

# Join predictions
percChangePreds <- rbind(colonySizePredictions, colonyNumberPredictions, areaPredictions)
percChangePreds$vary <- factor(percChangePreds$vary, levels = levels(d$vary))


percChangePreds$transmission <- factor(percChangePreds$transmission, labels = c('beta == 0.1', '0.2', '0.3'))


ggplot(d, aes(x = ValueChange, y = mean, colour = vary, shape = vary)) +
  geom_line(data = percChangePreds, aes(x = ValueChange, y = .fitted, colour = vary)) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = .04,
                position = position_dodge(width = 0.04)) +
  geom_point(size = 2.1, position = position_dodge(width = 0.04)) +
  ylab('Prob. Invasion') +
  xlab('Relative change') +
  scale_colour_poke(name = 'Population-level trait', 
                    labels = c( 'Group size', 'Number of groups', 'Host density'),
                    pokemon = 'oddish',
                    spread = 3) +
  scale_shape_manual(name = 'Population-level trait', 
                     labels = c( 'Group size', 'Number of groups', 'Host density'),
										 values = 15:17) +
  scale_y_continuous(breaks = c(0, 0.5, 1), 
                     labels = c('0.0', '0.5', '1.0'),
                     limits = c(0, 1)) +
  scale_x_continuous(breaks = c(0.25, 0.5, 1, 2, 4), 
                     labels = c(0.25, 0.5, 1, 2, 4), 
                     trans = 'log2') +
  facet_grid(transmission ~ ., labeller = label_parsed) +
	theme(legend.position = 'bottom', legend.direction = 'vertical')


%%end.rcode

\subsection*{Population size and host density}
The estimated GLM coefficients were always larger for simulations where population size was varied (sets \emph{ii} and \emph{iii}) than when host density (set \emph{i}) was varied (table~\ref{t:regrCoefs}, electronic supplementary material figure~S3).
Increasing population size, either by increasing group size or number of groups, increased the probability of invasion (figure~\ref{fig:plotValueChangeMeans}).
The positive relationship between group size and invasion was strong and significant at all transmission rates, while the relationship between number of groups and invasion was weaker but still significant.
In contrast, varying host density did not significantly alter invasion probability except for when $\beta = 0.3$ where there was a significant decrease in invasion probability with increased host density (GLM: coefficient = $\rinline{areaModel$estimate[6]}$, $p$ = \rinline{areaModel$p.value[6]}).
The 95\% confidence intervals for the coefficients of group size did not overlap with those for the coefficients of host density at any value of $\beta$ while the 95\% confidence intervals for coefficients of number of groups only overlapped with those for host density at $\beta = 0.2$.

\subsection*{Group size and number of groups}
In all cases, the estimated GLM coefficients were larger for simulations where group size was varied (set \emph{ii}) than when the number of groups (set \emph{iii}) was varied (table~\ref{t:regrCoefs}, electronic supplementary material figures~S3 and S4).
Increasing either group size or the number of groups increased the probability of invasion but this effect was  stronger for group size (figure~\ref{fig:plotValueChangeMeans}).
The 95\% confidence intervals of regression coefficients for group size did not overlap with those for number of groups for the simulations with $\beta = 0.2$ or $0.3$ but the confidence intervals did overlap for the simulations when $\beta = 0.1$.

%At the default parameter settings, the probability of invasion and persistence of the second pathogen was rare (figure~\ref{fig:plotValueChangeMeans} and electronic supplementary material, tables~S2--4).
%These proportions significantly increased with transmission rate (GLM: coefficient = \rinline{transTest$estimate[2]}, $p$ = \rinline{transTest$p.value[2]}).

\begin{table}
\captionsetup{width=0.9\textwidth}
\caption{Regression results}
\label{t:regrCoefs}
\centering
\begin{tabular}{@{}rlrrr@{}}
\toprule
$\beta$ & \emph{Variable} & \emph{Coefficient} & (95\% \emph{CI}) & $p$\\
\midrule
0.1   &    Group size   & \rinline{colonySizeModel$estimate[2]} & (\rinline{colonySizeModel$conf.low[2]}, \rinline{colonySizeModel$conf.high[2]}) & $< 10^{-4}$ \\
&    Number of groups & \rinline{sprintf('%.2f', colonyNumberModel$estimate[2])} & (\rinline{colonyNumberModel$conf.low[2]}, \rinline{colonyNumberModel$conf.high[2]}) & \rinline{colonyNumberModel$p.value[2]} \\
&    Density       & \rinline{sprintf('%.2f', areaModel$estimate[2])} & (\rinline{areaModel$conf.low[2]}, \rinline{sprintf('%.2f', areaModel$conf.high[2])}) & \rinline{sprintf('%.2f', areaModel$p.value[2])} \\[1em]
  0.2   &    Group size   & \rinline{colonySizeModel$estimate[4]} & (\rinline{colonySizeModel$conf.low[4]}, \rinline{colonySizeModel$conf.high[4]}) & $< 10^{-4}$ \\
&    Number of groups & \rinline{colonyNumberModel$estimate[4]} & (\rinline{colonyNumberModel$conf.low[4]}, \rinline{colonyNumberModel$conf.high[4]}) & \rinline{colonyNumberModel$p.value[4]} \\
&    Density       & \rinline{areaModel$estimate[4]} & (\rinline{areaModel$conf.low[4]}, \rinline{areaModel$conf.high[4]}) & \rinline{areaModel$p.value[4]} \\[1em]
0.3   &    Group size   & \rinline{colonySizeModel$estimate[6]} & (\rinline{colonySizeModel$conf.low[6]}, \rinline{colonySizeModel$conf.high[6]}) & $< 10^{-4}$ \\
&    Number of groups & \rinline{colonyNumberModel$estimate[6]} & (\rinline{colonyNumberModel$conf.low[6]}, \rinline{colonyNumberModel$conf.high[6]}) & $< 10^{-4}$  \\
&    Density       & \rinline{areaModel$estimate[6]} & (\rinline{areaModel$conf.low[6]}, \rinline{areaModel$conf.high[6]}) & \rinline{areaModel$p.value[6]} \\
\bottomrule


\end{tabular}
\caption*{
Regression results comparing effects of group size, number of groups and host density.
Estimated regression coefficients and their 95\% confidence intervals are from binomial GLM regressions with the proportion of invasions as the response variable and all explanatory variables being standardised by dividing by the default parameter value and applying a $\log_2$ transform.
Group size and number of groups were varied while keeping density equal while density was varied by changing geographic range size while keeping population size equal.
Results are given for three transmission ($\beta$) values.
}
\end{table}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section*{Discussion}

%todo what are the four things i want to say.
% Restate. Interpret at simple level.
% Comparative studies have problems and should consider these resiults
  % Further correlations between these variables.
% Limitations:
  % Two pathogens
  % Scope limits
% Global change affects these factors differently.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\tmpsection*{Restate results}

% Population-level traits such as increased increased host density \cite{morand1998density, kamiya2014determines, nunn2003comparative}, geographic range size \cite{nunn2003comparative, turmelle2009correlates, kamiya2014determines} and population structure (nonrandom interactions between individuals) \cite{maganga2014bat, turmelle2009correlates} have been shown to correlate with high pathogen richness, although the relationship with group sizes in many studies has been equivocal \cite{vitone2004body, gay2014parasite, rifkin2012animals, nunn2003comparative}. 

Overall, our results suggest that population size strongly promotes pathogen richness.
In contrast, host density was found to weakly promote pathogen richness if at all.
This is in agreement with theoretical arguments that population size is the more natural description of populations \cite{begon2002clarification} but is in contrast to the many comparative studies that find host density to be a significant predictor of pathogen richness \cite{morand1998density, kamiya2014determines, nunn2003comparative}.
This suggests that in these comparative studies, host density is acting as a proxy for population size rather than being a causal factor.
Alternatively, the local measurements of host density could be indicative of some other trait such as maximum host density, rather than global density.
Our results also suggest that a population made up of large groups would have higher pathogen richness than a population made up of many smaller groups.
Comparative studies examining group size have had conflicting results \cite{vitone2004body, gay2014parasite, nunn2003comparative, turmelle2009correlates} with a meta-analysis finding a non-significant relationship between group size and pathogen richness \cite{rifkin2012animals}.
The largest study specifically on bats found a negative relationship between group size and pathogen richness \cite{gay2014parasite}.
Therefore while our results regarding group size are potentially in conflict with the empirical literature, more definitive empirical studies are needed.
This conflict may also indicate a difference between pathogen richness accumulation (the focus of this study) and total pathogen richness.

Our results also suggest that host geographic range size does not promote pathogen richness, yet a number of studies have found evidence of this relationship \cite{kamiya2014determines, nunn2003comparative}.
In studies that do not also include population size or host density, geographic range size may be acting as a proxy for population size.
Alternatively this empirical association may be because in wild species, increased host geographic size tends to entail a greater variety of environmental conditions and a greater number of sympatric species which is known to also affect pathogen richness \cite{brierley2016quantifying} and these factors are not considered in this study.
Finally, we found an unexpected negative association between density and invasion probability at $\beta = 0.3$.
Due to the small size of the estimated coefficient, the marginal p-value and the lack of a consistent relationship at other $\beta$ values, we suspect that this is not a real effect but merely due to chance.
However, given that in our study, increased density affects disease dynamics by decreasing population structure, this negative relationship does fit with studies that suggest that population structure should allow pathogen coexistence \cite{maganga2014bat, turmelle2009correlates, allen2004sis}.

%The difference in size between these two effects might decrease in systems with higher dispersal rates.
%However, at very high dispersal rates, the population will be effectively unstructured and therefore changes in population structure will again have no affect on pathogen competition; it seems likely that population structure will be most important at intermediate dispersal rates.
%Similarly, the global host density of the species had little effect on pathogen invasion rate.
%In these simulations, increasing density without increasing population size was only achieved by reducing range size; this simply increased the number of connections between groups in the metapopulation network.
%This, in turn, increased the pool of susceptibles that were could be reached by one dispersal event of the invading pathogen.
%However, again, this effect was very weak compared to the strong changes in local disease dynamics caused by increasing colony size.


%\subsection*{Comparative studies}

Many comparative studies measure population-level traits, yet it is rarely considered how these might be causally related (though statistical correlations between explanatory variables are often handled appropriately).
For example, host density is often used in studies \cite{morand1998density, lindenfors2007parasite, nunn2003comparative}, yet host density is directly associated with population size.
Our results suggest that despite the association between these two traits, they are not equivalent.
%Our results suggest that it is in fact population size that is important, and previous authors have suggested that population size is the more natural measure to describe populations than host density \cite{begon2002clarification}.
%Therefore, density could be acting as a proxy for population size in these comparative studies.
These causal relationships between population-level traits should be considered more carefully in future comparative studies.
Researchers could include of an interaction term between geographic range size and host density to test for the importance of population size, for example.
%Furthermore, measurements of host density from the literature should be treated carefully; for example, were the original authors estimating global density or measuring local density in an area of unusually high host density?


% \tmpsection*{Only two pathogens, direct transmission}

This study was limited to one mechanism by which pathogen richness can be increased; the invasion and persistence of a newly evolved pathogen \cite{vickery1998parasite}.
However, other processes such as pathogen extinction are also likely to be important \cite{vickery1998parasite}.
We also restricted ourselves to the context of competition between two pathogens in a social host species.
As the model was of a directly transmitted pathogen there is no transmission between individuals in separate groups.
Infection via shared food sources or contact between individuals from different groups (e.g., during swarming \cite{kerth2008causes}) would act to reduce population structure and therefore host density might become more important. 
However, further modelling would be required to demonstrate this while only empirical studies would be able to indicate the true relative importance of these different transmission routes in wild populations.


%\subsection*{Global change}

It is clear that many species are suffering strong population changes due to global change \cite{thomas2004extinction}.
These changes might affect geographic range size \cite{thomas2004extinction}, population size \cite{craigie2010large}, population connectivity \cite{rivera2015habitat, fonturbel2014forest} or group size \cite{lehmann2010apes, atwood2006influence} to different degrees.
The monitoring of these different aspects of population change, especially in bats, can often be difficult and may require further developments in monitoring to be effective, for example developing methods that use data from acoustic detectors \cite{walters2013challenges, measey2016counting, lucas2015generalised}.
Our results suggest that pathogen communities will respond differently depending on which traits are most strongly affected by global change.
In short, species suffering reductions in groups size \cite{lehmann2010apes, atwood2006influence} are predicted to experience a decrease in pathogen richness in the long term and there is some evidence that this process is occurring \cite{altizer2007threatened, turmelle2009correlates}.
Species that are experiencing an increase in group size \cite{lehmann2010apes} would be expected to gain new pathogen species and therefore pose a greater risk of being the source of a zoonotic disease.
In contrast, species suffering a reduction in geographic range size \cite{thomas2004extinction} or a decrease in population size \cite{craigie2010large}, without a corresponding decrease in group size,  are expected to experience smaller changes in pathogen richness.
Only by examining the mechanisms that control pathogen richness can we understand and predict these changes.

\section*{Data accessibility}
\sloppy
The implementation of the model is available as an R package on GitHub \cite{metapopepi}. This can be found at \href{https://github.com/timcdlucas/MetapopEpi}{https://github.com/timcdlucas/MetapopEpi}.
All other code and simulation output data is available on GitHub at \href{https://github.com/timcdlucas/Abundance-Density-Manuscript}{https://github.com/timcdlucas/Abundance-Density-Manuscript}.
\fussy


\section*{Competing interests}
We have no competing interests.



\section*{Author's contributions}

TCDL wrote the simulations and performed the analysis.
TCDL, HMH and KEJ co-designed the study.
TCDL drafted the manuscript.
TCDL, HMH and KEJ all edited the manuscript and gave final approval for publication.

\section*{Acknowledgements}
We thank Andy Fenton and David Murrell for comments on the manuscript and help with the analytical model. 
This work, The Dynamic Drivers of Disease in Africa Consortium, NERC project no. NE-J001570-1 was funded with support from the Ecosystem Services for Poverty Alleviation Programme (ESPA). 
The ESPA programme is funded by the Department for International Development (DFID), the Economic and Social Research Council (ESRC) and the Natural Environment Research Council (NERC).


\section*{Funding}
This study was funded through a CoMPLEX PhD studentship at University College London supported by BBSRC and EPSRC (TCDL). KEJ was funded by the Ecosystem Services for Poverty Alleviation Programme (ESPA) (NE-J001570-1).





\bibliography{epilit}
\bibliographystyle{amnat}

\include{abundance-density-tl-SI}



\end{document}