Models with environmental drivers offer a plausible mechanism for the rapid spread of infectious disease outbreaks in marine organisms

The first signs of sea star wasting disease (SSWD) epidemic occurred in just few months in 2013 along the entire North American Pacific coast. Disease dynamics did not manifest as the typical travelling wave of reaction-diffusion epidemiological model, suggesting that other environmental factors might have played some role. To help explore how external factors might trigger disease, we built a coupled oceanographic-epidemiological model and contrasted three hypotheses on the influence of temperature on disease transmission and pathogenicity. Models that linked mortality to sea surface temperature gave patterns more consistent with observed data on sea star wasting disease, which suggests that environmental stress could explain why some marine diseases seem to spread so fast and have region-wide impacts on host populations.

where Ki is carrying capacity within the cell and Ri is a constant daily recruitment. Both cellspecific values are set to default values (K and R), then scaled by region (Table C2). with vj,t indicating mean along-shore current near cell j at time t, cell width cw, and standard deviation θQ. The dispersal kernel used here combines a long-range asymmetric transport process (the along-shore current) with a short-range Gaussian diffusion-advection process (the specific destination cell post-transport). Note that Eq. B3 evaluates all possible source cells for a destination cell and thus Eq. B4 inverts the typical dispersal calculation. To simplify, we assume that the propagule source is a single point at the midpoint of the source cell and that the current at the source cell represents conditions for the entire time step. We then integrate across the destination cell width to convert distance into cell-to-cell transport. Our model makes the assumption that a virionor other infectious agentis not absorbed after it contacts a host (see below for implications).
Transition from E to I occurs within a cell at a constant proportion pBG (=0 for the purely environmentally-driven model) or when cumulative degree-days DDi,t approach or exceed a threshold: The scaling exponent c governs the degree to which partial transition occurs as stress levels approach the threshold. Degree-days accumulate exponentially above a local threshold for extreme temperatures: where a and b are scaling parameters, Ai,t indicates temperature anomaly in cell i at time t, Tmin indicates the minimum absolute temperature for temperature anomalies to be considered stressful, and DDrecover indicates the maximum amount of accumulated stress recovered on nonanomalous days. Alternative forms of the stress function are detailed in Appendix D.
Both Ti,t and Ai,t, the daily mean sea surface temperature and temperature anomaly in cell i at time t, are calculated using ROMS data for the Pacific coast from January 1st, 2013 through December 31st, 2015 49 . Anomaly is defined as deviation from local mean seasonal cycle as determined using 1999-2011 climate data. We linearized the data as described in Appendix A and averaged across gaps to match the spatial layout of the model (2500 km coastline in 5 km cells). We used the same approach to transform the daily along-shore surface current vi,t from the ROMS model.

Const-Q-SEI model
The Const-Q-SEI model is a special case of the EnvDr-Q-SEI model for which pBG=1.0 and τ=20. Consequently, exposed S class individuals develop a full infection after 20 days, transforming to E class, then transition to I after a single timestep.

EnvInf-Q-SEI model
The where ∆Ti,t is a 5-time-step moving variance average scaled such that the mean daily value across the simulation =0.5, or 50% of the transmission rate for Const-Q-SEI, but ranges from close to 0 on cold days to 10x or greater on the hottest days.
Secondly, disease progression is accelerated during periods of high temperature 39,41 . For Const-Q-SEI, there is a constant delay transitioning from S to E which, for the default value of τ=20, always takes 20 timesteps to complete. For EnvInf-Q-SEI, however, higher temperatures speed up this transition: disease progression doubles for every 3 °C above the Tcold threshold, to a maximum of 4x faster at 6+ °C. Thus, a single high temperature day advances the disease 2-4x as fast as a cold day, and the complete S->E transition could take anywhere from 5-20 timesteps depending on the proportion of hot and cold days in the incubation period following initial infection.

Virion loss
Virion loss upon sea star contact is likely a relatively small source of pathogen loss compared to other sources. Furthermore, because each cell in our model has similar host density, there would be no variation in relative loss rate due to this mechanism. If we were modeling filter feeders, loss due to host contact might be quite high, causing us to over-estimate infection rates. This assumption would have been most problematic for dose-dependent exposure in a filter feeder with variable density 66 . In such a case, the assumed loss of virions upon infection in dense host populations could keep virion densities below the infective dose, preventing disease spread.
Such inverse density dependence could set up spatially patchy mass mortalities that impact low density host populations, but this scenario does not appear to fit SSWD.

Appendix D: Sensitivity to stress, outbreak, and model parameters
Note: in all referenced parameter sensitivity tables except

Stress-accumulation parameters
In our model, we make two broad assumptions about how temperature-related stress accumulation affects sea stars. Firstly, as shown in Equation B.6, we assume that stress accumulates non-linearly using an exponential function (repeated here for convenience): where a and b are scaling parameters. Although some studies of intertidal organisms have shown a non-linear response to temperature (e.g., oyster mortality 57 and mussel growth 58 ), the precise form of the relationship is unclear and may differ for sea stars. Coral bleaching, for example, is often modeled with stress as constant above a threshold (i.e., a step-function). Consequently, we tested several different possible stress functions to assess model sensitivity to the mechanism of stress accumulation ( Figure D1). We found that the overall goodness of fit for the EnvDr-Q-SEI was relatively unaffected by the shape of the stress function (Table D1).
Secondly, we assumed that accumulated stress only affected the transition from the asymptomatic E class to the symptomatic I class. In reality, stress is likely to increase mortality for the uninfected S class as well. To test the importance of this assumption, we altered the constant survival terms σS and σE as follows: where x indicates class, m is a constant stress-related mortality modifier, and DDi,t is the accumulated degree-days within cell i at time t. By varying m, we showed that increasing stressrelated mortality shifted the model with minimal MLE values from EnvInf-Q-SEI to EnvDr-Q-SEI and, at very high levels, to Const-Q-SEI (Table D2).

Parameter estimation
We used MLE values to select model parameters that were hard to estimate from field data or experimental work. We chose the infection parameter β (Table D3), infection delay τ (Table D4), and initial infection location (Table D5) based on which value best balanced likelihoods for the density and prevalence data. Note that relative model support, and thus our overall results, stays consistent across most parameter values. In Table D6 we compare model likelihoods at different data aggregation levels, and show that the models are closest at the medium resolution we chose. The EnvInf-Q-SEI model is favored in both categories at finer aggregation and the EnvDr-Q-SEI model has a better prevalence MLE value at coarser aggregation. Note that MLEs cannot be compared across different aggregation levels.

Hypercube sampling of model parameters
Because this is an exploratory rather than explanatory model, we did not use specific real-world disease parameters and were instead focused on a qualitative comparison between the output and the SSWD data. However, the results could still depend on the specific parameters chosen rather than the structure of the model. We estimated demographic parameters for P.
ochraceus from Menge 67 and set epidemiological parameters from experimental results and by selecting parameters to minimize the MLE values (Tables D3-5). To assess sensitivity to changes in the demographic and epidemiological parameters of the core disease model, we ran a Latin hypercube sampling analysis 68 , allowing parameters to freely vary by +/-25% (or across a specified range; see Table C1). The results from this analysis show that, although variance is introduced and there is a wider spread for the possible MLE values (especially for the density MLE and Const-Q-SEI; Figure D2a), the relative support for each model variant remains consistent across almost all parameter combinations ( Figure D2b) and our overall conclusions remain the same. Note the relative lack of points in the lower-left corner of Fig. D2b (the "negative/negative" region), which shows that there were very few instances for which Const-Q-SEI was better on both metrics than the more environmentally-dependent models (though the reverse is often true), or EnvInf-Q-SEI than EnvDr-Q-SEI. In most instances, there was a tradeoff between density and prevalence likelihood.
We used partial rank correlation coefficient (PRCC) analysis to determine which parameters significantly affected the relative difference between the models (Table D7)   The x-axis indicates the difference in density MLE values between the two models indicated, and the y-axis indicates the difference in prevalence MLE values. A positive value means that the first model was less supported by the data than the second model.