Table 1 provides a summary of the six non-human primate studies conducted at US Army Medical Research Institute of Infectious Diseases (USAMRIID). All animal studies were conducted under an Institutional Animal Care and Use Committee (IACUC) approved protocol in compliance with the Animal Welfare Act, PHS Policy, and other Federal statutes regulations relating to animals and experiments involving animals; the research facility is accredited by the International Association for Assessment and Accreditation of Laboratory Animal Care and adheres to principles stated in the Guide for the Care and Use of Laboratory Animals, and the National Research Council, 2011.

In studies 1 and 2, rhesus and cynomolgus macaques, respectively, were exposed to the Marburg Angola virus referred to as MARV (Marburg virus/H.sapiens-tc/ANG/2005/Angola-1379c – USAMRIID challenge stock “R17214”). In study 3, cynomolgus macaques were exposed to Ebola virus (EBOV) at a target dose of 100 plaque forming units (pfu; Ebola virus/H.sapiens-tc/COD/1995/Kikwit-9510621; 7U EBOV; USAMRIID challenge stock “R4415;” GenBank # KT762962). In study 4, African green monkeys were exposed to the Malaysian Strain of Nipah virus (NiV) isolated from a patient from the 1998 to 1999 outbreak in Malaysia and provided to USAMRIID by the Centers for Disease Control and Prevention. In study 5, cynomolgus macaques were exposed to the Josiah strain of the Lassa virus (LASV; challenge stock “AIMS 17294;” GenBank #s JN650517.1, JN650518.1). In study 6, African green monkeys were exposed to Yersenia pestis (Y. pestis), a causative bacterial agent of bubonic and pneumonic plague. Additional details on studies 1, 5, and 6 have been published elsewhere (Malhotra et al., 2013; Johnston et al., 2015; Ewers et al., 2016). Dependent on the study, animals were exposed under sedation via either aerosol, intramuscular (IM) injection, or intratracheal (IT) exposure (see Table 1).

In each study, the animals were implanted with remote telemetry devices (Konigsberg Instruments, Inc., T27F or T37F, or Data Sciences International Inc. L11: see details in Table 1) 3–5 months before exposure, and, if used, with a central venous catheter 2–4 weeks before. They were then transferred into BSL-3 (bacterial exposures) or BSL-4 (viral exposures) containment 5–7 days prior to challenge. Baseline data from the telemetry devices were collected for 3–7 days before exposure. Monitoring via the telemetry devices continued until death or the completion of the study. The mean duration of recordings within each study is shown in Table 1. The exposure time (t = 0) denotes the time of IM injection or IT exposure or when a subject was returned to the cage following aerosol exposure (~20 min).

Data from a total of 46 animal subjects from the six studies were available. Eight subjects from the EBOV exposure study were excluded from post-exposure analysis because they received therapeutic interventions following the challenge, which could be a confounding factor. Five subjects across the cohorts were excluded on the basis of either substantial data loss from equipment failure or development of fever more than 2 days prior to the studies’ mean (i.e., possible co-morbid infections or complications). This resulted in a total of 33 animal subjects for our analytical cohort (N = 33).

Multimodal physiological data from the animal subjects were made available in NSS format (Notocord Systems, Croissy-sur-Seine, France). The multimodal physiological data from the implanted telemetry devices included raw waveforms of the aortic blood pressure (sampling frequency f_s = 250 Hz), electrocardiogram (ECG; f_s = 500 Hz), intrathoracic pressure (f_s = 250 Hz), and core temperature (f_s = 50 Hz). All signals were measured internal to the animals, which generally resulted in very high-signal fidelity. Using Notocord software, we extracted the features listed in Table 2 from the raw waveforms.

We categorized all features retrieved from the Notocord system (Table 2) as: pre-exposure baseline or simply baseline – data collected from the start of the recording up to 12 h before the viral or bacterial challenge; and post-exposure – data collected 24 h after the viral or bacterial challenge until death or the completion of the study. Relative timing of the labeled regions is depicted in Figure 2A. Data from 12 h before and 24 h after viral or bacterial challenge were excluded from performance metrics due to differences in animal handling and exposure sedation that resulted in significant physiological deviations from baseline data unrelated to pathogen infection.

All subjects across all six studies developed fever as a result of the pathogen exposure. For early warning, fever onset is an important reference point. We define fever onset as the first time that the subject’s core temperature measurement exceeds 1.5°C above that subject’s diurnal baseline (Laupland, 2009) with the additional constraint that the temperature is sustained above threshold for at least 2 h. Leveraging this fever onset, we further categorized the post-exposure data as being pre-fever – before onset of fever; or post-fever – after onset of fever.

Figure 2B illustrates the key steps taken to condition the features (Table 2) extracted from the Notocord system after the data labeling. These steps are applied per subject and per feature to reduce diurnal and inter-subject time dependencies in the data. As an example, the time series feature, Bazett, which represents QT-corrected intervals from an ECG waveform are shown at each processing stage.

We train our random forest models on two post-exposure stages, thus allowing the algorithms to adapt to different physiological cues during the pre-fever and post-fever phases. The pre-fever random forest model is optimized to discriminate the earliest stages of illness by training it on pre-fever data samples vs. baseline. The post-fever random forest learns discriminants of the febrile phase of illness by training it on post-fever data samples vs. baseline. The number of data points used for training is balanced for equal representation of the classes.

Both models, pre-fever and post-fever, are trained using the l-minute epoch summary statistics generated in the data pre-processing step, including mean, 25^th and 75^th quantiles for all 12 features listed in Table 2.

The models are implemented using the TreeBagger class in the MATLAB Statistics and Machine Learning Toolbox.

We next apply a two-stage detection process, depicted in Figure 2C, to the prediction scores generated by the pre- and post-fever models for each l-minute epoch, with the primary goals of reducing the overall false alarm rate and incorporating recent historical scores in the decision.

In stage one of the detection process, a time series of feature vectors is processed on two parallel paths. One path calculates a pre-fever random forest score while the other path independently calculates a post-fever random forest score. On each path, the score is compared to a threshold associated with the respective model (threshold selection is described in section “Model Tuning”). Initial detections occur when a score exceeds the threshold. To reduce the likelihood of spurious detections, we buffer the initial detections over a window of n epochs and perform binary integration (Shnidman, 1998), calculating a moving average over an l*n minute window.

In stage 2 of the detection process, the parallel paths are reunited by taking the maximum of the pre-fever and post-fever moving average value at each epoch. This combined score is compared to a second stage threshold of m/n, where m is an integer such that m ≤ n. Combined scores in excess of m/n are declared to be in the exposed class, and we use the term “declaration” to denote the final decision from the two-stage processing. Note that the buffering aspect of binary integration imposes some latency on the system, so no declarations are reported in the first l*n minutes.

We evaluate overall performance of our models using three key performance metrics: probability of detection, P_d, probability of false alarm, P_fa, area under the receiver operating curve (AUC), and mean early warning time Δt.

We calculate the probability of correct declaration, P_d, as the number of true positive declarations over the total number of post-exposure samples. In addition, we compute P_d on the subset of pre-fever and post-fever samples. We use the term, system P_d, to represent correct detection over all post-exposure data samples regardless of fever status, while pre-fever P_d indicates the refinement, where correct detections are evaluated exclusively on the subset of pre-fever data samples. The probability of false alarm, P_fa, (also referred to as the system P_fa) is defined as the number of false positive declarations over the total number of baseline samples. In order to estimate small false alarm rates with meaningful precision, we require a large number of baseline data samples. For false alarm analysis, we supplement with baseline data from some animals that were excluded from the primary analysis. These data include seven full days of measurements from each of nine animals prior to pathogen exposure: seven subjects from the EBOV study (excluded due to therapeutic intervention following exposure) and two subjects from the NiV study (which developed fever earlier than our exclusion criteria). We compute 95% confidence intervals for P_d and P_fa assuming normal distribution as the number of trials is large (the number varies depending on the dataset and metric under evaluation, but is greater than 500 for all scenarios considered here).

We generate receiver operating characteristic (ROC) curves to measure system performance by calculating P_d vs. P_fa at a series of threshold values (sweeping the first-stage detection threshold while holding the second-stage m/n threshold constant) and report the AUCs evaluated against pre-fever and post-fever data samples.

Another important measure of system performance is the mean early warning time, Δt. The early warning time for an individual subject is defined as the time of the first true declaration minus the time of fever onset. We compute the mean across all subjects to characterize the early warning time afforded by the system and report 95% confidence interval based on a t-distribution since the number of subjects is small (<30 for the two subgroups considered here).

We evaluate detection performance under three distinct scenarios to address our core research questions. First, to answer the fundamental question of how well pathogen exposure can be detected based solely on physiological measurements, we focus on data from the subset of N = 20 animal subjects from the EBOV and MARV studies (Studies 1–3). We develop our algorithms using a 3-fold cross-validation approach, which has been shown to perform better (Shao, 1993) than leave-one-out validations for small dataset. This approach explicitly varies five experimental variables (species and sex of animal, exposure route, pathogen, and target dose; see Table 1) across the three partitions, which reduces the likelihood of biasing the model for any particular condition.

Second, to evaluate whether the early warning capabilities of our algorithm extend to other pathogens, we train the models on the N = 20 subjects from EBOV and MARV and apply them to an independent dataset of N = 13 animal subjects from the LASV, NiV, and Y. pestis studies (Studies 4–6).

For the third research question, we investigate performance degradation when the inputs of the classifier are restricted to emulate a limited set of measurements that could be obtained with a non-invasive wearable device. While direct measurement of aortic blood pressure, core temperature, and intrathoracic pressure rely on invasive or intrusive sensors, ECG signals can be readily measured with wearable sensors. For this evaluation, we limit the classifier inputs to the set of EGG-derived features. This scenario is also trained on the EBOV and MARV data set and applied to the LASV, NiV and Y. pestis data set.

Model tuning, including feature selection and other classifier and detection parameters is performed using systematic parameter sweeps within the subset of N = 20 animals exposed to EBOV and MARV. Table 3 summarizes the tunable parameters from sections Data Pre-Processing: Conditioning Physiological Data, Random Forest Ensemble, and Detection Logic.

To characterize performance as a function of different epoch lengths and number of trees, we make use of the random forest out-of-bag errors. Random forest ensembles are generated through an aggregated bagging process, whereby a random subset or “bag” of data points is selected with replacement to build a decision tree. The process is repeated until a specified number of trees are generated. Out-of-bag errors are calculated during training by evaluating decision trees against the samples that were not in their bag, providing a convenient assessment of classifier performance. We sweep the parameter values for k, l, and n_trees and consider tradeoffs for both the out-of-bag errors and the computation times for pre-processing the data.

For feature selection and determining the number of features, n_features, we assign one of the three cross-validation partitions for parameter tuning and the remaining partitions for model training and performance validation. We use a backward elimination feature selection method, leveraging the out-of-bag errors to iteratively identify and drop the feature ranked least important. The impact of varying n_feature is further characterized by the P_fa and pre-fever P_d in this parameter tuning partition.

First-stage detection thresholds are selected based on a user-defined target P_fa, evaluated at the final stage of the declaration logic. We sweep the first-stage thresholds independently for the pre-fever and post-fever classifier and select the smallest threshold for each model that achieves a target P_fa ≤ 0.01.

For the second-stage detection parameters, we sweep m and n, from n = 1 (30 min) to n = 36 (18 h) and m = 1, 2,…, n. For each pair (m, n), we evaluate performance in terms of the mean early warning time Δt and AUC. We also consider the performance when m is set to the estimated optimal threshold for a constant (non-fluctuating) signal in noise, mopt≈10−0.02n0.8 (Shnidman, 1998).

We begin by evaluating trade-offs for the k- and l-minute epoch length. Figure 3A shows the relative computation time for the data conditioning steps of section “Data Pre-Processing: Conditioning Physiological Data” as a function of the normalization epoch, k-minutes. Preprocessing was performed on a Dell desktop computer with dual Intel Xeon E5607 processors and 12GB RAM. Preprocessing is very time-consuming for the shortest epochs but the time burden decreases with increasing epoch length, leveling out around 30 min. Figure 3B shows the impact of both the feature normalization and summary statistics epoch lengths on classifier performance. In general, shorter epochs for the feature normalization are associated with lower errors and therefore better classification accuracy. In contrast, long epochs for the summary statistics provide better accuracy than short ones. The result suggests selecting l ≥ k, that is, the summary statistics epoch should be at least as long if not longer than the normalization interval. As a balance between processing time and classifier accuracy, we select our epochs as k = l = 30 min. With 48 epochs in a 24 h period and three summary statistics per physiological features listed in, we nominally compute 144 features per subject per day.

Next, we optimize the random forest parameters, n_trees and n_features. As shown in Figure 4, both the false positive rate and the pre-fever Pd improve as n_features increases from 1 to about 10, but performance plateaus beyond 10. Similarly, classifier accuracy improves as n_trees increases but plateaus beyond about 15. We settle on a classifier composed of 15 trees grown on the 10 highest ranked features.

Table 4 shows the 10 highest-ranked features for each of the three partitions used in cross validation as well as the 10 highest-ranked features when we re-train the models on the full set of N = 20 animals. The ranked feature importance shows consistency with clinical symptomology, namely that core temperature-based features (mean and quantiles of temperature) in the post-fever model rank highest in importance. Before fever, however, ECG- and blood pressure-derived features are among the highest in feature importance, as has been reported at the earliest stages of sepsis (Korach et al., 2001; Chen and Kuo, 2007; Ahmad et al., 2009; Scheff et al., 2012). Among the blood pressure features, quantiles of systolic and diastolic aortic pressure rank as the most important. Among ECG-derived features, means and quantiles of QT intervals [corrected (Bazett, 1920; Fridericia, 2003) or not], RR intervals and PR intervals are routinely selected as those with the greatest predictive capability. Respiratory rate features derived from the intrathoracic pressure waveform were seldom ranked among the most important.

The first stage threshold values, listed in Table 3, are estimated for the pre- and post-fever models to enforce a target P_fa = 0.01. Using these first-stage thresholds, we see in Figure 5, that larger n (longer buffer length) enables slightly better detection capability in the sense of AUC, but at the expense of reduced early warning time. The estimated optimal threshold, denoted by the dashed line, is reasonably aligned with peak performance for both early warning time and AUC at each n, allowing for a methodical assignment of the second-stage threshold, m_opt, given a binary integration window, n. In this analysis, we select n = 24 and m = m_opt = 11 to achieve high AUC while maintaining a system latency of no more than 12 h.

A time series of the combined score resulting from the two-stage detection process for a representative animal subject from the MARV aerosol study is shown in Figure 6A. The combined score, for this subject, remains below the detection threshold (dashed horizontal line) before virus challenge, rises sharply around exposure (which is excluded) due to anesthesia, then rises again at ~2 days post-exposure, where the first “exposed” declaration (dashed vertical green line) occurs.

In this cross-validation assessment, we evaluated performance over a total of 9,931 decision points from N = 20 subjects and found a system P_d = 0.80 ± 0.01, a pre-fever P_d = 0.56 ± 0.02, a system P_fa = 0.01 ± 0.003, and a mean early warning time of Δt_mean = 51 ± 12 h. Detailed performance metrics on each subject can be found in the Supplemental Material.

We further evaluated algorithm performance for all subjects with the family of ROC curves shown in Figure 6B, where the P_d is separately evaluated against pre- and post-fever data samples. For this three-fold cross-validation, we find AUC = 0.93 ± 0.01 for pre-fever data, and AUC = 0.99 ± 0.001 for post-fever data.

Figure 6C shows a plot of correct declarations as a function of early warning time. This plot focuses on detectability in the pre-fever region for a threshold corresponding to P_fa = 0.01. Mean early warning time, estimated for each pathogen exposure is shown as a dashed vertical line, which indicates individual differences between pathogens and exposure study conditions. Among these three studies, we see the earliest mean warning time for MARV IM exposure at Δt_mean = 69 ± 16 h, while the two aerosol exposures, EBOV and MARV, have similar mean values at Δt_mean = 33 ± 26 h and Δt_mean = 39 ± 18 h, respectively.

We test our pre- and post-fever models against data from the LASV aerosol, NiV intratracheal, and Y. pestis aerosol studies (Table 1, N = 13 subjects). The combined score vs. time is shown in Figure 7 for one representative subject for each pathogen.

This independent validation set includes over 11,000 decision points including the supplemental baseline data from nine subjects that were otherwise excluded (the supplemental points contribute only toward P_fa; they are excluded from P_d calculations). The corresponding ROC curves and mean early warning times for the independent validation set are shown in Figure 8A. Again, detection performance against post-fever samples is nearly perfect, and we observe significant pre-fever positive predictive value of the model, with an AUC = 0.95 ± 0.01. Across the three pathogens, we find a system P_d = 0.90 ± 0.007 and P_fa = 0.03 ± 0.004, a pre-fever P_d = 0.55 ± 0.03, and a mean early warning time of Δt_mean = 51 ± 14 h. Calculating Δt_mean for each pathogen exposure study individually, we find that the NiV IT study has the longest Δt_mean = 75 ± 30 h (though NiV subjects also have the longest incubation period, ~5 days), and that LASV aerosol and Y. pestis aerosol exposure studies have Δt_mean = 33 ± 26 h and Δt_mean = 41 ± 25 h, respectively (with a mean incubation period ~3.5 days). A summary of the performance metrics from this independent validation data set are shown along with the cross-validation data set performance in Table 5.

As an in silico simulation for degrading our animal dataset to what may be collected using a wearable monitoring device for humans, we reduced the considered feature set to include only ECG-derived features such as RR, QT, QRS, and PR intervals. Figure 8 compares our algorithm performance using all available features (Figure 8A) and features derived only from the ECG waveform (Figure 8B). For the degraded feature set, we see only modest decreases in performance with Δt_mean = 46 ± 14 h, pre-fever P_d = 0.55 ± 0.03, and system P_d = 0.89 ± 0.008 and P_fa = 0.03 ± 0.004.