Abstract
Postinjury epilepsy is an potentially preventable sequela in as many as 20% of patients with brain insults. For these cases biomarkers of epileptogenesis are critical to facilitate identification of patients at high-risk of developing epilepsy and to introduce effective anti-epileptogenic interventions. Here, we demonstrate that delayed brain–heart coincidences serve as a reliable biomarker. In a murine model of post-infection acquired epilepsy, we used long-term simultaneous measurements of the brain activity via electroencephalography and autonomic cardiac activity via electrocardiography, in male mice, to quantitatively track brain–heart interactions during epileptogenesis. We find that abnormal cortical discharges precede abnormal fluctuations in the cardiac rhythm at the resolution of single beat-to-beat intervals. The delayed brain–heart coincidence is detectable as early as the onset of chronic measurements, 2–14 weeks before the first seizure, only in animals that become epileptic, and increases during epileptogenesis. Therefore, delayed brain–heart coincidence serves as a biomarker of epileptogenesis and could be used for phenotyping, diagnostic, and therapeutic purposes.
SIGNIFICANCE STATEMENT No biomarker that readily predicts and tracks epileptogenesis currently exists for the wide range of human acquired epilepsies. Here, we used long-term measurements of brain and heart activity in a mouse model of post-infection acquired epilepsy to investigate the potential of brain–heart interaction as a biomarker of epileptogenesis. We found that delayed coincidences from brain to heart can clearly separate the mice that became epileptic from those that did not weeks before development of epilepsy. Our findings allow for phenotyping and tracking of epileptogenesis in this and likely other models of acquired epilepsy. Such capability is critical for efficient adjunctive treatment development and for tracking the efficacy of such treatments.
Introduction
No biomarker that readily predicts and tracks epileptogenesis (Gowers, 1881) currently exists for the wide range of human acquired epilepsies, developed after brain insults such as traumatic brain injuries, stroke, and infections (Engel et al., 2013; Pitkänen et al., 2016a,b). These biomarkers are critical to identify patients at risk of epilepsy, introduce early and effective interventions to prevent establishment of epilepsy, and evaluate therapeutic efficacy of such treatments.
Many investigators have focused on chronic cardiac alterations and autonomic dysfunction as biomarkers for epileptogenesis (Kheiri et al., 2012) as well as features for seizure prediction (Massé et al., 2013; Kolsal et al., 2014; Moridani and Farhadi, 2017; Pavei et al., 2017) with limited success. However, the underlying coupled brain–heart dynamics and the alterations in those dynamics during epileptogenesis are not fully investigated.
According to the lockstep phenomenon (LSP; Lathers et al., 1987), epileptiform discharges can induce intermittent synchronized cardiac sympathetic and vagal neural discharges. LSP is therefore a representation of seizure-induced deficits in the brain circuitry of the cardiovascular system (Goodman et al., 2008) that may alter, although not fully impair, the system's function. We hypothesized that these events are the extreme end of smaller, potentially abnormal, fluctuations in brain activity that are likely to subtly modify cardiac function long before the microscale instability propagates and results into an overwhelming seizure (Bragin et al., 2000; Schevon et al., 2008; Stead et al., 2010). Our aim was therefore to investigate the frequency and strength of the signals from brain affecting heart as a function of progression of epileptogenesis. Once detected and quantified the coupled brain–heart dynamic could potentially serve as a prognostic and diagnostic biomarker of epileptogenesis.
Previously, we developed a murine model of post-cerebral malaria epilepsy that duplicates elements of the human post-infection acquired epilepsy. As with the human condition, the animal model presents with long and highly variable epileptogenesis periods (Ssentongo et al., 2017) thus provides a platform to investigate the potential of brain–heart interaction as a relevant biomarker for epileptogenesis.
Here, we continuously track brain–heart dynamics during epileptogenesis. To achieve this, we developed a new metric of cardiac rhythmicity whose statistics are independent of behavior and state of vigilance (SOV). We report that in animals that survive the malarial infection and become epileptic, extrema of the cardiac metric are preceded by potentially abnormal fluctuations in brain activity. The causal transmission of activity along the brain–heart axis starts early and progresses during epileptogenesis. These findings can thus be used to identify potential patients at risk of developing epilepsy, track the progression of epileptogenesis as well as anti-epileptogenic treatments, and predict the occurrence of seizures.
Materials and Methods
All animal work was approved by and performed under administration of the Institutional Animal Care and Use Committee at the Pennsylvania State University.
Experimental design
This work uses data acquired in the development of a murine model of post-cerebral malaria (CM) epilepsy previously described in (Ssentongo et al., 2017). We examined combinations of mouse and parasite strains for development of post-CM epilepsy (Ssentongo et al., 2017): male Swiss Webster, male C57BL/6, and male CBA/CaJ mice and Plasmodium berghei ANKA (PbANKA) and Plasmodium berghei NK65 (PbNK65) parasites. Detailed methods are included in the study by Ssentongo et al. (2017). Critical elements for this work follow.
Surgical procedure and care.
Mice were implanted with hippocampal, cortical, and electromyography (EMG) electrodes at least 5 d post-treatment following the procedures described by Ssentongo et al. (2017). A cohort of the animals also received electrocardiography (ECG) electrodes to monitor cardiac activity. The details of the electrode fabrication and implant are discussed by Ssentongo et al. (2017). At the completion of the surgery, animals were returned to their individual home cages for recovery. Mice were allowed to recover and monitored for a minimum of 3 d post-surgery. Animals were then cabled and recordings with continuous periods of at least 2 weeks lasted until the animals expired.
Data analysis
At least 1 week of data was analyzed for each animal. Longer datasets up to when the animals developed their first seizure, lost the ECG lead, or were euthanized were used when available. Our inclusion criteria for each hour of data involved availability of viable ECG and at least one low-noise hippocampal depth and electrocorticography (ECoG) recordings.
All recorded data were inspected via in-house written LabVIEW (National Instruments) and MATLAB (MathWorks) programs that allow for simultaneous re-referencing, filtering, spectral analysis and annotation.
Line length.
Line length (LL) for discrete time series is calculated from 〈|xn − xn−1|〉τw, where τw denotes the averaging window length. The averaging is achieved with a low-pass infinite impulse response (IIR) filter with cutoff frequency of 1/τw, applied both forward and backward through the data to eliminate phase delays.
Cortical discharges.
Cortical discharges were detected by a custom script implemented in MATLAB (MathWorks). The raw electroencephalography (EEG) data were first bandpass filtered using 15–250 Hz pass band to preserve the high-frequency nature of the discharges while eliminating the low-frequency artifacts. LL of the EEG (EEGLL) was then computed using a window length τw = 300 ms. Cortical discharges vary in pattern, amplitude, and duration throughout epileptogenesis. The window length was selected to capture potentially abnormal EEG events as well as the variety of epileptic discharges including sharp waves, spikes, and spike and wave complexes (Kane et al., 2017).
State of vigilance.
State of vigilance (SOV) was scored using a semiautomatic approach via linear discriminant analysis (LDA) similar to the algorithm described by Sunderam et al. (2007). Briefly, for each animal, 4–6 h of video-EEG data within 1 d were manually scored for SOV. A training set was then generated with EEG spectral power in frequency bands 0.5–4, 4–8, 8–12, 12–25, and 25–80 Hz, plus neck EMG power used as features. The remaining 18–20 h of data were set aside as out-of-sample test data. We marked onset of transitions between different SOVs: rapid eye movement (REM), non-rapid eye movement (NREM), and wake.
We further verified the performance of the classifier for the out-of-sample test data and used the classified SOV as the new training set to recursively update the parameters of the LDA. For most animals under analysis the overall automated accuracy exceeded 90%.
Heart rate characterization
QRS complex.
In each heartbeat, the QRS complex was detected via a custom written MATLAB script (MathWorks). The ECG is bandpass filtered from 2 to 250 Hz to remove the low (DC drift) and high-frequency artifacts. LL with window size τw = 50 ms was used to highlight the QRS complexes (∼40 ms long) from the background ECG activity. The standard deviation (SD) of the ECGLL time series was estimated in non-overlapping 2-s-long windows. The peak point of the QRS complex, the R wave position, was then detected via threshold crossing of 2–3 times the SD from the mean as the corresponding location of the local maximum. The point-by-point RR interval (RRI) and heart rate time series are then calculated based on the detected R wave times.
Motion artifacts and other glitches not removed with filtering are detected by convolving the unfiltered ECG with a step function. Any QRS complex detected within 20 ms of those glitches is then marked as a potential false detection and is excluded from analysis. We removed hour-long blocks of data with >2 glitches from the analysis. For each detected glitch one QRS complex and 2 RRIs are excluded. Therefore, for the hour-long blocks that are included in the analysis, the maximum number of excluded RRIs is 4 of 18,000 RRIs, which is negligible.
Cardiac metric dlnRR.
The cardiac metric dlnRR was calculated at any time, T, over window lengths 2τ long:
where
where fractional number of intervals include partial intervals that remain inside the window. Therefore, the discretization effects are minimized.
In the limit as τ → 0, if the derivative of RR time series is defined, this metric becomes the derivative of the logarithm of the RRI:
dlnRR is a local measure that is proportional to the normalized RRI in the window length over which it is calculated. Therefore, the window length τ directly affects dlnRR's capability to efficiently track and represent the underlying RRI time series. If τ is too small, then dlnRR is dominated by normal fluctuations between consecutive RRIs induced by motion or cardio–respiratory coupling. If τ is too large, then the numerator becomes insensitive to single RRI fluctuations, while the denominator is too long to follow state-dependent changes in the RRI time series. However, in Figure 3D we show that choice of τ does not affect the state-independence of the dlnRR.
Statistical analysis
Our analysis seeks to quantify correlations between fluctuations in the extrema of brain dynamics, expressed in EEGLL, and fluctuations in the extrema of cardiac activity, from dlnRR. The distributions of EEGLL and dlnRR are not Gaussian; therefore, we used nonparametric rank-order statistics. The ranks for each metric were computed and normalized over either 1-h- or 1-d-long blocks of data approximating a uniform selection of ranks on the range of (0,1). As detailed in Results, we use the Wilcoxon rank sum test for samples of either metric conditioned on values of the other. Under an assumption of stochastic independence of these metrics, the normalized rank (rank sum/N) of such a sample of size N (for N ≫ 100) should have mean = 0.5, and SD = 
. All rank-ordered statistics are then presented either in absolute units (see Fig. 4) or as mean and SD of the rank-sum distributions (see Figs. 5, 6). Expected false detection rates from null hypotheses with p < 0.005 are identified.
Results
We examined a combination of mouse and parasite strains for development of post-cerebral malaria epilepsy (Ssentongo et al., 2017). Animals were continuously recorded from after ample postsurgical recovery time for at least 2 weeks or until they expired.
From the animals that survived to recordings, 75% developed epilepsy as defined by observation of two or more seizures with durations >10 s and clear behavioral manifestations (Ssentongo et al., 2017). The epileptic animals all had relatively long latencies to the first seizure with median of 5–10 weeks postinfection (Fig. 1A) depending on mouse/parasite combination. Further, we observed incidents of sudden unexplained death in epilepsy (SUDEP) and other seizure-related mortalities such as gradual decline of physiological signs, heart rate and behavior, after multiple seizures in 10% of the epileptic cohort.
Chronic recording and analysis summary. A, Listed are the distribution of animals chronically recorded that met the criteria for analysis of brain–heart coincidences during epileptogenesis. Numbers in parentheses indicate total number of mice studied for development of epilepsy post cerebral malaria. For detailed description of model statistics and mouse–parasite combinations see Ssentongo et al. (2017). These animals showed similar seizure rates and latencies to the larger cohorts studied for development of the post-cerebral malaria epilepsy model (Ssentongo et al., 2017). B, Animals were implanted and continuously monitored (orange lines) after treatment and ample recovery time. For animals that became epileptic, the latency to the first convulsive seizure (marked with magenta asterisks) ranged from 39 to 115 d post-infection (shaded gray area). Because the physiological measurements did not start until the animals fully recovered from cerebral malaria, we assume that we only captured part of the epileptogenesis process. Instances of substantially long single RRIs (Subs. LRRI) were present in all epileptic mice. The onset of appearance of such events is marked by black open circles. The black dashed lines separate epileptic, non-epileptic, and control mice. Non-epileptic mice were rescued from cerebral malaria and chronically recorded from but did not develop seizures. Data before the first convulsive seizure for epileptic mice (shaded gray area), and until ECG leads broke in control (entire recording time) and non-epileptic mice (dark green crosses) were analyzed. Data from 30 and 73 d post-infection from animal #9 are used in later Figures 4 and 5.
To investigate the brain–heart coupling during epileptogenesis, nine epileptic mice with viable ECG electrodes, at least until their first seizure, were selected. This subset is representative of the original epileptic cohort both in terms of seizure rate and distribution as well as latency to first convulsive seizure (Fig. 1A), which ranged from 39 to 115 d postinoculation (Fig. 1B). In addition to the nine epileptic mice, four animals that were rescued from cerebral malaria but did not become epileptic (termed non-epileptic mice), and four control mice (treated but not infected) were investigated. Animals were excluded from analysis when they developed spontaneous recurrent seizures or their ECG lead became unviable.
Single substantially long RRIs are preceded by epileptic discharges
Cardiac arrhythmia manifested in visually distinguishable, substantially long RRIs in the period between implant and the occurrence of the first convulsive seizure. RRI is defined as the interval between successive R waves marked by the peak of the QRS complex in each heartbeat. We quantitatively identified substantially long RRIs by a user-defined criterion of any interval >150 ms, which is ∼150% of the mean RR value in mice. We found that high-frequency cortical and hippocampal discharges precede instances of long RRI by at least 500 ms. An example of such precedence is shown in Figure 2A.
Substantially long RRIs are preceded by epileptiform cortical discharges. Instances of abnormal cortical and hippocampal discharges preceding substantially long RRIs are shown for one animal. Substantially long RRIs are detected via the user-defined threshold of any interval >150 ms. A, High amplitude, high-frequency EEG activity immediately precedes the long RRI (red trace in ECG) onset of which is marked by the gray vertical bar. Time-series traces of hippocampal activity are recorded from left and right dorsal hippocampi [depth hippocampus left (DL,) and depth hippocampus right (DR)], ECoG recorded from frontal and S1 cortices [ECoG left anterior (ELA), ECoG right anterior (ERA), ECoG left posterior (ELP), ECoG right posterior (ERP)]. ECG is bandpass-filtered between 15 and 250 Hz to eliminate DC artifacts and overall trend of the signal. B, In the same event as in A, the LL of the brain activity at the right posterior cortical site (LLRP), computed over 300-ms-long windows, increases preceding a long RR interval. C, Almost all substantially long RRIs (31 of 35) are preceded by large increases in the LL of the right posterior ECoG shown in color-coded values. Overlaid is the average of the LL of the right posterior ECoG activity (blue solid line) for all detected abnormally long RRIs along with the 95% confidence bounds (dashed blue lines). The data extracted from animal #3 (Fig. 1B), A and B, are from 50 d post-infection, 11:00 A.M.
LL of the EEG (EEGLL), calculated with window lengths τw = 300 ms, for the right posterior cortical electrode increases immediately before the rise in the RRI time series (Fig. 2B) for the event shown in Figure 2A. Similar increases in EEGLL were observed in the rest of the cortical and hippocampal channels. For the animal presented in Figure 2 this particular temporal pattern was observed for the majority (31 of 35) of substantially long RRIs (Fig. 2C). Inspection of the associated recording video confirmed that these events were subconvulsive.
A limited number of substantially long RRIs was observed for each animal that became epileptic. This ranged from 10 to 53 instances of long RRIs of the 7–84 million RRIs marked from the recording times before the animals' first seizures. Almost all of these cases were preceded by subconvulsive cortical and hippocampal discharges reflected in increased EEGLL.
The substantially long RRIs, >150 ms, were observed starting many weeks into the epileptogenesis process (Fig. 1B, black open circles). This led us to question whether earlier in epileptogenesis abnormal brain activity would induce more subtle changes in duration of single RRIs and whether the size of such an effect would continuously increase before the first seizure. Detection of abnormal transmission early in epileptogenesis requires (1) a cardiac measure that identifies small abnormal fluctuations in the RRI, and (2) the statistical power to correlate these RRI fluctuations with extremely small fluctuations in brain activity that otherwise would not be identified as epileptic.
State-independent cardiac metric: dlnRR
Distributions of conventional metrics of cardiac activity such as absolute duration of RRIs or heart rate variability (HRV) measures vary with behavioral states and are highly sensitive to transitions between such states (Baharav et al., 1995; Ronkainen et al., 2005; Boudreau et al., 2013; Hajek and Buchanan, 2016). Detrending methods such as detrended fluctuation analysis (Peng et al., 1995) attempt to eliminate the behavioral state dependency. These methods however fail to quantify nonstationarities at the resolution of single RRIs (Bryce and Sprague, 2012; Govindan et al., 2014). HRV measurements conditioned on one State of Vigilance (SOV) are also prone to error due to the challenges in identifying quiet, resting states and transitions in and out of them. Therefore, to investigate the relationship between brain activity and cardiac rhythm, we developed a new and more flexible cardiac metric sensitive to small but abnormal fluctuations in the duration of single RRIs.
We defined the cardiac metric dlnRR(T) at time T proportional to the difference in mean RRI before and after T, computed in windows τ wide, normalized by the average RRI over the period 2τ. We denote this metric dlnRR(T) because it is the discrete estimate of the logarithmic derivative of the RRI.
The distribution of dlnRR is approximately symmetric with near-zero mean. Unlike conventional measures of cardiac activity such as heart rate or HRV parameters (i.e., root mean square of successive differences (RMSSD) of the RRIs), distribution of dlnRR is not influenced by changes in behavioral state (Fig. 3A,B).
Temporal measures of cardiac rhythm for a control mouse. A, Distributions of dlnRR, RMSSD, and RRI for different SOV pooled from 2 weeks of recordings from control animal #2. Distributions of dlnRR and RMSSD are calculated using identical window lengths (τ) of 0.5 s. RMSSD and RRIs have varying distributions across different states of vigilance; two-sample Kolmogorov–Smirnov test, REM and NREM (pRR = 6.41–6, pRMSSD = 9.1e-4), REM and wake (pRR = 5.9e−8, pRMSSD = 5.95e−4), and NREM and Wake (pRR = 3.9e−4, pRMSSD = 1.12e−4). The distribution of dlnRR remains similar and approximately symmetric with near-zero mean for all states of vigilance. By construction, for the τ = 0.5 s, single long RRIs lead to values in the upper tail of the dlnRR distribution. Dashed gray line indicates the boundary for the 99% upper tail of the dlnRR distribution. B, Mean and SD of RRIs, RMSSD, and dlnRR distributions across different SOVs. C, Mean (solid lines) and SD (dashed lines) of RMSSD calculated over different window lengths (τ). For all τ, within the physiologically relevant range of 5–20 RRIs, mean and SD of RMSSD distributions for different SOV are different. D, Mean (solid lines) and SD (dashed lines) of dlnRR calculated over a sweep of τ remain relatively similar for all SOV.
Sensitivity of dlnRR to choice of τ is discussed in the Materials and Methods section. For our analyses we selected τ = 0.5 s to calculate dlnRR, mean heart rate, and RMSSD. For τ = 0.5 single long RRIs that occur within the window after time T result in large positive values of dlnRR in the tail of the distribution.
Delayed coincidences between brain and heart fluctuations
To investigate correlations between changes in RRIs and small fluctuations in brain activity, in each 1 h block of data, LL of the EEG and dlnRR for each beat time were calculated. LL of one hippocampal lead was then selected for further analysis of the brain–heart coupling. We expect that if abnormal cortical activity is transmitted from brain to heart, then we would find coincidences with positive time delay between potentially abnormal values of EEGLL and potentially abnormal values of dlnRR.
To identify transmission along brain–heart axis, our analysis leverages statistical tests of dependence of the EEGLL and dlnRR distributions. By construction, potentially abnormal events will have values in the tails of each of these distributions. Note that by choosing simply the tails for the EEGLL, for most of the epileptogenesis period, the events detected are too small and ill-defined in shape to identify as clearly abnormal or epileptic. In addition, by construction, the number of potentially abnormal events per hour is constant, as defined by the detection threshold applied, and independent of time since implant or of the animal's health.
The joint distributions of EEGLL and dlnRR ranks, P(EEGLL(t), dlnRR (t + Δt)), over different time offsets, Δt, for an hour-long block of data are shown in Figure 4A. Except for Δt = 1 s, these joint distributions are consistent with independence of the individual metrics. At Δt = 1 s, the joint distribution is highly non-uniform with high densities in the extrema of the ranks as indicated by darker corners.
Joint rank distributions of brain (EEGLL) and cardiac metric (dlnRR) over different time offsets (Δt) for an hour-long block of data. The metrics are first transformed into the rank space. A, The joint EEGLL-dlnRR rank distributions are then calculated for different time offsets (Δt). When Δt = 1 s, joint rank distribution of EEGLL and dlnRR is no longer uniform such that higher ranks of dlnRR correspond to higher ranks of EEGLL. The normalized CRS of EEGLL (left) and dlnRR (bottom) highlight the non-uniformity of the joint distribution at Δt = 1 s. The CRSs of EEGLL and dlnRR (black line) increase for higher ranks and decrease for lower ones. In contrast, the mean rank-sum for P (rank EEGLL (t), rank dlnRR (t+Δt)|Δt=10 s) remains at ∼0.5 (red line) which indicates a uniform distribution. Dashed gray lines indicate p value = 0.1. B, Shown are the CRSs of EEGLL as a function of Δt for the same animal as in A at different times post-infection (PI). Contour lines mark select iso-probability (yellow = 0.1, black = 0.2) of the Wilcoxon rank sum test with the null hypothesis that the ranks are chosen from independent distributions. Outside the time range 0 ≤ Δt ≤ 2 s, this test has PWilcoxon ≫ 0.2, which is consistent with independent brain and heart metrics. In contrast, within the small range of 0 ≤ Δt ≤ 2 s the CRS diverges from independence for the extrema of the ranks (magenta and cyan areas). That p ≤ 0.1 (A, gray dashed lines) inside the 0 ≤ Δt ≤ 2 s interval indicates a dependence between large fluctuations in brain or cardiac activity with a time delay of 0–2 s. The brain–heart interdependence is smaller early (30 d PI) and worsens during epileptogenesis (73 d PI). dlnRR and EEGLL extracted from epileptic animal #9 (Fig. 1), 30 and 73 d PI, 12:00–1:00 P.M.
We quantified the non-uniformity of these joint distributions by computing the conditional rank-sum (CRS) of EEGLL and dlnRR rank distributions. These conditional distributions are normalized to have expected mean rank of 0.5. The CRSs are shown for offset times Δt = 1 s (black) and Δt = 10 s (red), for the rank-sum of dlnRR as a function of EEGLL (Fig. 4A, bottom inset) and for the rank-sum of EEGLL as a function of dlnRR (Fig. 4A, left inset). For independence of EEGLL and dlnRR we expect a uniform distribution and CRS = 0.5, as observed for Δt = 10 s. In contrast, the CRS at Δt = 1 s diverges away from that, consistent with interdependent distributions.
To confirm the dependence of the joint distribution on the time offset (Δt), we calculated the CRS of EEGLL and dlnRR rank distributions over a sweep of time-offsets (−50 ≤ Δt ≤ 50 s). In mice that later became epileptic the CRS of EEGLL (Fig. 4B) and dlnRR (data not shown) during epileptogenesis diverged significantly from uniformity (CRS = 0.5) within the time range of 0 ≤ Δt ≤ 2 s. Significance here is quantified from the Wilcoxon rank sum test p value, which is the probability that the observed rank-sums would have come from random rank samples. Further, in animals that later became epileptic, the statistical interdependence of the EEGLL and dlnRR distributions progressively strengthened during epileptogenesis. This pattern is indicated in smaller CRSs within 0 ≤ Δt ≤ 2 s earlier during epileptogenesis (Fig. 4B, bottom; 30 d post-infection) and larger CRSs within 0 ≤ Δt ≤ 2 s later during epileptogenesis (Fig. 4B, top; 73 d post-infection). Overlaid on Figure 4B are the iso-contours for p value of 0.1 (yellow) and 0.2 (black). We will use repeated independent excursions through this contour to yield highly significant observations (p < 0.005).
Quantification with the Spearman's rank-order correlation yields similar confirmation of statistical interdependence of EEGLL and dlnRR distributions that peaks with a time-delay of ∼1 s (data not shown).
Brain–heart coincidence as a biomarker for epileptogenesis
The significant divergence of the joint dlnRR-EEGLL distribution from independence with a peak delay of 1 s implies that fluctuations in brain activity potentially cause or lead to fluctuations in cardiac activity. Because we observed this trend and its progression only in mice that later became epileptic, we hypothesized that the significant divergence of the joint dlnRR and EEGLL distribution from uniformity within the specific range of 0 ≤ Δt ≤ 2 s could be used to distinguish animals that would become epileptic from the ones that did not. These brain–heart coincidences, once quantified, would serve as a biomarker of the epileptogenesis.
The apparent brain to heart coupling is primarily mediated by isolated abnormal events. By construction these events fall in the upper tail of the EEGLL and dlnRR distributions. To further quantify these events during epileptogenesis, we define a conditional rank threshold (cRT) and select EEGLL and dlnRR ranks that are above this threshold. The cRT = 0.99 is marked by the gray dashed line in Figure 3A for the dlnRR distribution. We then extract values of the other metric as a function of time with respect to these high-rank values, brain or heart triggered events, and quantify, as a function of time offset, whether their distribution is consistent with independence. This is achieved by computing the p value of the Wilcoxon rank sum test, which we denote as Psig.
Optimization of Psig to classify epileptic and non-epileptic mice
We adopted a cross-validation classification method to optimize Psig to separate epileptic from non-epileptic and control animals. We first computed the CRS for EEGLL (with cRT = 99%) for every hour-long block of data for all animals. We then selected every other day of measures from the first 14 d of recordings for all animals as a training set, and found the minimum CRS that maximized the classification specificity for identifying epileptic animals. The Wilcoxon p value associated with this CRS is Psig. The CRS threshold was then applied to all remaining data, and hours with CRSs larger than this threshold were marked for existence of brain–heart coincidences.
The fraction of hours per day with significant brain–heart coincidences with the optimized CRS threshold (Psig = 0.0554) is shown in Figure 5A. Because each day's value is constituted from 24 independent measures, the probability of observing a fraction of hours ≥0.25 with the null hypothesis of random distribution of ranks, is <0.002. For the epileptic animals, these fractions fluctuate during the first 21 d after recording onset, but they all generally increase during epileptogenesis. In contrast, for non-epileptic animals they remain low from the onset of recording until time-of-death.
Evolution of brain–heart coupling in epileptic, non-epileptic, and control mice using hour-long blocks of data. A, Daily fraction of hours marked for existence of brain–heart coincidence are shown for epileptic (colored diamonds), non-epileptic (colored squares), and control (hexagons) mice; based on a cRT of 0.99. The fractions increase during epileptogenesis with large fluctuations within the first 21 d of the recordings. In contrast, they remain low for animals that were rescued from cerebral malaria but did not become epileptic (non-epileptic) and controls. The hours are marked for existence of brain–heart coincidences if they have conditional EEGLL rank sums larger than a specified rank-sum threshold. The rank-sum threshold is selected to maximize classification specificity for a collection of every other days from the first 14 d of recordings. Mice were then classified according to a classification boundary on fraction of marked hours to maximize specificity (black solid line). This boundary is represented by the gray dashed line separating epileptic, non-epileptic, and control distributions in B (top). B, Bottom, The ROC curve using the boundary and fraction of marked hours. All epileptic mice but one (A, light blue trace) are separated from non-epileptic and control mice after 21 d of recording. C, Specificity (blue circles), sensitivity (orange crosses), and the minimum conditional rank sums (gray triangles) selected as the classification boundary to mark hours for brain–heart coincidences, as a function of cRT. End point for epileptic animals is the occurrence of the first spontaneous seizure and for non-epileptic and control cohorts is time of death.
We then defined a classification criterion by maximizing specificity (Fig. 5A, black line, B, gray dashed line), using the fraction of marked hours per day. Under this criterion we can correctly classify all 13 mice under analysis into epileptic (n = 9) and non-epileptic (n = 4) groups over almost all of the period shown. Only one epileptic animal is temporarily misclassified within the third week of recording.
Optimization of cRT
We further investigated the sensitivity of the rank-sum classifier to our criteria of potentially abnormal events, cRT. Shown in Figure 5C are the minimum CRS thresholds selected as the classification boundary to mark hours for brain–heart coincidences over a sweep of conditional rank thresholds [cRT ∈ (85%, 99%)]. As more ranks are allowed into the potentially abnormal samples, the CRS threshold (and Psig) increases to separate epileptic from non-epileptic and control animals with maximum specificity. Maximum specificity is achieved at the expense of degrading sensitivity. But, that neither specificity nor sensitivity decreases below 80% indicates the classification is robust over a wide range of cRTs.
Day-long analysis blocks
By expanding our analysis block size from 1 h blocks to full-day intervals we increased the statistical power of the analysis. As done previously, we computed CRS of EEGLL with cRT = 0.99 for full day-long blocks of data, and then derived the maximum CRS within the window of 0 ≤ Δt ≤ 2 s.
The maximum CRS for EEGLL within the range 0 ≤ Δt ≤ 2 s is shown in Figure 6A for each day for all epileptic, non-epileptic, and control animals. As with the daily fractional detections (Fig. 5A), this metric steadily increases during epileptogenesis (for the epileptic animals) with smaller fluctuations during the first weeks of recordings. With this statistic the epileptic animals are completely separated from non-epileptic animals from the first day of recording (Fig. 6B). Therefore, the metric provides a classifier with complete sensitivity and specificity.
Evolution of brain–heart coupling in epileptic, non-epileptic, and control mice using day-long blocks of data. A, CRS of the cardiac-triggered brain events from each day for all animals: epileptic mice (colored diamonds), non-epileptic mice (colored squares), and control mice (hexagons); based on cRT of 0.99. The color-coding for all mice is the same as in Figure 5A. Similar to daily fractional detections in Figure 5A, this statistic increases during epileptogenesis for epileptic mice while it remains low for non-epileptic and control mice. All epileptic animals, including the animal that was transiently misclassified in Figure 5A, are fully separable from non-epileptic and control mice as shown in the distributions of these animals in the top of B (green, epileptic; black, non-epileptic; magenta, controls). Therefore, this statistic provides a classifier with complete sensitivity and specificity (shown via the ROC curve in B, bottom). Epileptic Animal 9 in Figure 4 is marked by the purple trace. A, Inset, The number of epileptic (gray), non-epileptic (green), and control (magenta) alive at days post-recording onset. C, Specificity (blue circles) and sensitivity (orange crosses), as a function of cRT. The sensitivity values are computed based on a classification boundary to maximize specificity. As the cRT is relaxed to include more intermediary values of the rank[dlnRR(t)] distribution, the sensitivity and specificity decrease and the classifier performs worse. Gray triangles indicate the minimum classification threshold that maximizes specificity in the units of fraction of SD away from the mean of the control and non-epileptic distributions. End point for epileptic animals is the occurrence of the first spontaneous seizure and for non-epileptic and control cohorts is time of death.
We further note that each day's measurement constitutes an independent measure, with one free variable; cRT. This analysis therefore separates 480 independent measures from animals that became epileptic, from 308 independent measures including ones from non-epileptic animals (168) and controls (140).
The potentially abnormal events in the analysis in Figure 6A are based on a cRT = 0.99. To evaluate the sensitivity of the classifier presented in Figure 6B to this threshold, we repeated the analyses for a sweep of cRTs over (75–99%). The animals were then separated with a classification boundary to maximize specificity. As shown in Figure 6C, the classification sensitivity begins to degrade for cRTs <95%, although specificity decreases only when the cRT is <88%. Again, the analysis is quite robust for a broad range of statistically relevant thresholds over the rank distributions.
Qualitatively similar results can be achieved without the inherent threshold introduced by detection of potentially abnormal events through cRT, by using Spearman's rank correlation analysis (data not shown).
Univariate measures do not track epileptogenesis
As shown in Figure 7A, the mean EEGLL remained relatively stationary in the duration before the first seizure in epileptic mice and during entire lifetime for non-epileptic mice. Similarly the SD of dlnRR distribution although slowly increased for a portion of epileptic animals (n = 4), did not differentiate epileptic from non-epileptic mice.
dlnRR and EEGLL alone do not separate epileptic from non-epileptic cohorts. Distribution of mean EEGLL (A) or SD of dlnRR (B) as a function of time is not different for epileptic (colored squares) versus non-epileptic (brown hexagons) animals. Time (t − τsz) is with respect to the first seizure for the epileptic mice and death for non-epileptic mice. Even hours before the first seizure (right) the mean values are not differentiable between epileptic and non-epileptic mice. The distributions of mean of EEGLL (C) and SD of dlnRR (D) are pooled from 10 d intervals of recordings from all epileptic (magenta) and non-epileptic mice (green). The two distributions are overlapping which indicates that as a feature for classification neither of them provides sufficient sensitivity or specificity. ROC curves are calculated for a variety of statistical measures of the EEGLL (C, right) and dlnRR (D, right) distributions including mean, SD, SD normalized by the mean, and values of the 1–5% upper bounds. As shown by the ROCs, none provide significant sensitivity or specificity to separate the epileptic from non-epileptic animals.
Both EEGLL and dlnRR distributions significantly overlap for epileptic and non-epileptic animals. We tested a variety of statistical measures of each of the distributions of EEGLL and dlnRR as features to separate epileptic from non-epileptic and control mice. None provided significant sensitivity or specificity to separate, as illustrated in Figure 7, C and D, by the receiver operating characteristics (ROC) curves of mean, SD, SEM, and values of the 1–5% upper bounds.
Discussion
We used continuous long-term recordings in a murine model of post-cerebral malaria epilepsy to derive metrics of brain–heart interactions. We then investigated the statistical codependence of these metrics as a biomarker of the epileptogenesis period in nine mice that later developed epilepsy under the acquired epilepsy model. Underlying these metrics are fluctuations in brain activity causing, ∼1 s later, fluctuations in cardiac rhythmicity. We used LL as the metric for brain activity to specifically highlight epileptiform-like events and transient cortical discharges. We developed a SOV-independent metric of cardiac rhythmicity to highlight abrupt changes in cardiac rhythm with resolution of single long RRIs. Transitions between SOVs induce transient changes in RRIs. The dlnRR is therefore designed to be state-independent to not confound the state induced changes with abnormal variations of RRIs. By construction, abnormally long brain and cardiac events result in large positive values of the EEGLL and dlnRR.
We demonstrate that abnormal cortical discharges precede abnormally long single RRIs. These delayed coincidences were only detected in animals that survived the malarial infection and later became epileptic. Further, once quantified, we observe that the strength of the coupling between abnormal brain and heart events increases during epileptogenesis.
Our observations that subconvulsive cortical discharges lead to isolated long RRIs implicate involvement of the autonomic nervous system impinging on heart. The mechanistic basis for this phenomenon might be imbalanced parasympathetic activity transmitted to the heart via the vagal nerve or a decrease in sympathetic activity. The detected positive time-delay of 1 s is consistent with previous reports of an induced reduction in heart rate after electrical stimulation of the vagal nerve branch innervating the heart (Rosenblueth and Simeone, 1934; Warner and Cox, 1964).
Impaired cardiac function is considered a contributing clinical manifestation of severe childhood malaria (Mockenhaupt et al., 2004; Yacoub et al., 2010; Herr et al., 2011). One might question whether our observations are due to post-malarial cardiac dysfunction. Animals that were identified as non-epileptic in this study experienced and then were rescued from cerebral malaria and therefore were subject to potential cardiac damage. Compared with animals that later became epileptic we found negligible incidents of abnormally long RRIs and brain-to-heart delayed coincidences in non-epileptic mice. Control animals presented with even lower incidence of development of abnormal transmission along brain–heart axis. Our findings are thus not attributed to the cardio-pathological effects of malarial infection.
In the context of epilepsy, our findings support the body of evidence that epileptic discharges and seizure activity recruit cortical structures as well as autonomic nuclei regulating cardiac rhythmicity (Schraeder and Lathers, 1983; Lathers et al., 1987; Oppenheimer and Cechetto, 1990; Oppenheimer et al., 1991). The coalescence models in acquired epilepsies argue that during epileptogenesis the propagation of abnormal activity across small populations forms larger interconnected epileptogenic networks (Bragin et al., 2000; Bikson et al., 2003; Schevon et al., 2008; Stead et al., 2010; Wang et al., 2017). We found that in animals that later become epileptic the brain–heart coincidence becomes stronger during epileptogenesis. This likely represents recruitment of additional autonomic structures into the pathologic circuitry originally impaired by the initial cerebral malaria insult. If we assume that this coalescence is the substrate for the epileptogenesis, then interventions to reverse it would lead to decreases in the brain–heart delayed coincidences.
The strength of the brain-to-heart coupling during epileptogenesis shown in Figures 5 and 6 exhibited an increasing pattern with day-by-day fluctuations in the first 3–4 weeks of recordings. Compared with the large fluctuations in the daily fractional detections (Fig. 5A), the day long detection metric (Fig. 6A) shows relatively smaller day-by-day variations. Therefore, the day-long detections separate pre-epileptic and non-epileptic animals with perfect sensitivity and specificity from the first day of the recordings. The high degree of separability of pre-epileptic from non-epileptic and control mice provided can thus be used in phenotyping animals in studies and assessing risk of epilepsy in patients.
The large fluctuations of the daily fractional metric (Fig. 5A) in all pre-epileptic mice reflect the finer temporal variations of the brain-to-heart coupling. These fluctuations were the most severe in the animal that was transiently misclassified (Fig. 5A, light blue trace). We associate these changes with the dynamic nature of epileptogenesis. We hypothesize that successful treatments to interrupt epileptogenesis will be reflected as decreases in the abnormal brain-to-heart coupling captured via the daily fractional detections.
We note that our measurements, although continuous and long-term, started weeks after the malarial infection. Therefore, the true onset of detectable neurophysiological changes during epileptogenesis remains unclear.
We assert that in future work, our observation of development of abnormal brain-to-heart coupling before the first convulsive seizure can be adapted to forecast seizure clusters with long seizure-free intervals between them. These clusters would be reflected in increases of the abnormal brain-to-heart delayed coincidences followed by gradual decreases representative of the seizure-free periods.
The cortically-induced disturbances of cardiac rhythm have the potential to make the network more susceptible to more seizures and vulnerable to their effects (Lathers et al., 1987; Devinsky et al., 1997; Hilz et al., 2002; Altenmüller et al., 2004; Schuele et al., 2007). The LSP is evidence that epileptic discharges can interrupt autonomic regulation. Therefore LSP is proposed to share a common underlying mechanism with SUDEP (Schraeder and Lathers, 1983; Leutmezer et al., 2003; Altenmüller et al., 2004; Dütsch et al., 2006; van der Lende et al., 2016). Our observation of abnormal brain-to-heart coupling implies that cortical discharges are involved in induced autonomic imbalance that could lead to catastrophic conditions such as gradual deterioration of physiological signs and increased risk of SUDEP. For these cases, the strength of the brain-to-heart coupling should be investigated to assess the need and efficacy of cardiac pacemakers to reduce risk of SUDEP.
Although there are quite a few animal models of epilepsy, the time course of disease development and the underlying processes do not necessarily match human epileptogenesis periods. The murine model of post-CM epilepsy studied here models the human conditions of post-infection acquired epilepsy. Both express long and variable epileptogenesis periods before observations of spontaneous unprovoked seizures. Our findings of brain–heart delayed coincidences in this model can offer a novel approach to prevent epilepsy in high-risk patients—post-traumatic brain injuries, post-infection, post-anoxic/ischemic, postsurgical—through pharmacological trials with biomarker monitoring.
A critical need in biomarker research is development of preferably noninvasive and less costly markers that can be applied easily to human population. We detected the brain-to-heart coupling during epileptogenesis from delayed coincidences between EEG and ECG. Electrocardiograms are non-invasive. But further investigation of the underlying mechanisms of brain–heart axis may require measurements from the vagal nerve innervating the heart. Direct measurements of the cardiac branch of the vagal nerve can be provided by bipolar electrodes and implanted radio transmitters (Jung et al., 2006; Ogawa et al., 2007; Tan et al., 2008; Choi et al., 2010; Sevcencu et al., 2016).
The brain–heart coincidence measure introduced here is based on fast temporal correlations, of the order of 1–2 s, between variations in single RRIs and small fluctuations in brain activity. The temporal resolution needed for detection of such correlation is enabled because mouse heart beats many times per second. The analyses presented here need to be adjusted for the relatively slower heart rates in human subjects for clinical applications.
The work presented here highlights one of the many potential mechanistic coincidences between CNS and much more accessible physiological measurements. These physiological measures such as cardiac activity for parasympathetic and skin sympathetic nerve activity for sympathetic nervous system (Doytchinova et al., 2017; Uradu et al., 2017), limb motion (Zafeiriou et al., 1999; Hellwig et al., 2000, 2001; Kanemaru et al., 2014), cerebrovascular function (Bar-Klein et al., 2017), and SOV (Sedigh-Sarvestani et al., 2014) cannot only provide us with more insight into the underlying mechanisms of neurologically-sourced pathologies but can lead to novel diagnostic approaches and identification of new phenotypic features.
Footnotes
This work was supported by Pennsylvania State Institute for the Neuroscience from Pennsylvania Department of Health Tobacco Funds, a Multidisciplinary Grant from Citizens United for Research in Epilepsy (CURE), National Institute of Health Grant R01EB019804, and a doctoral Academic Computing Fellowship from Pennsylvania State University to F.B. We thank Ali Nabi, Balaji Shanmugasundaram, and Myles W. Billard for assistance in development of the acquisition electronics and recording electrodes.
Conflict-of-interest: F.B., B.J.G., and S.J.S. are co-inventors on a U.S. provisional patent No. 62/663,213. The remaining author declares no competing financial interests.
- Correspondence should be addressed to either Fatemeh Bahari or Bruce J. Gluckman, Pennsylvania State University, University Park, PA 16802, fzb116{at}psu.edu or BruceGluckman{at}psu.edu
