## Abstract

Electrical communication between cardiomyocytes can be perturbed during arrhythmia, but these perturbations are not captured by conventional electrocardiographic metrics. We developed a theoretical framework to quantify electrical communication using information theory metrics in two-dimensional cell lattice models of cardiac excitation propagation. The time series generated by each cell was coarse-grained to 1 when excited or 0 when resting. The Shannon entropy for each cell was calculated from the time series during four clinically important heart rhythms: normal heartbeat, anatomical reentry, spiral reentry and multiple reentry. We also used mutual information to perform spatial profiling of communication during these cardiac arrhythmias. We found that information sharing between cells was spatially heterogeneous. In addition, cardiac arrhythmia significantly impacted information sharing within the heart. Entropy localized the path of the drifting core of spiral reentry, which could be an optimal target of therapeutic ablation. We conclude that information theory metrics can quantitatively assess electrical communication among cardiomyocytes. The traditional concept of the heart as a functional syncytium sharing electrical information cannot predict altered entropy and information sharing during complex arrhythmia. Information theory metrics may find clinical application in the identification of rhythm-specific treatments which are currently unmet by traditional electrocardiographic techniques.

## 1. Introduction

The human heart consists of 5 billion autonomous cardiomyocytes [1] with simple rules of operation and minimal central control. The behaviours of individual cardiomyocytes are orchestrated by electrical conduction between adjacent cells connected by specialized cell-to-cell junctions called intercalated discs [2]. The intercalated discs contain gap junctions with large non-selective connexin channels that allow ions and other small molecules to diffuse freely in the cytosol of adjacent cells and reduce internal electrical resistance [3]. By providing low-resistance connections between cardiomyocytes, gap junction channels allow cardiac excitation to propagate rapidly throughout the heart [4].

However, this task can be perturbed during abnormal heart rhythm (*cardiac arrhythmia*). Cardiac excitation is an electrical wave spatially spanning from a region of action potential depolarization (*wavefront*) to a region of action potential repolarization (*wavetail*) [5]. When a wavefront meets a wavetail, the electrical wave breaks at the intersection, and the electrical information flow can be interrupted at the wavebreak [6]. Furthermore, a completely different output can be generated by arrhythmic circuits where the wave rotates around the wavebreak. Therefore, cardiac arrhythmia can lead to loss of information about the input. Conventional electrocardiographic metrics can measure the sequence of electrical excitation [7–9], but cannot quantify how arrhythmia impacts the communication between individual cardiomyocytes.

Information theory is a mathematical theory of communication to quantify information [10]. Information theory has been successfully used to evaluate biological communication in computational neuroscience [11], transcriptional regulation [12,13], bacterial quorum sensing [14], chemotaxis [15], biochemical signalling networks [16] and evolutionary biology [17]. Information theory metrics such as mutual information can quantify the sharing of information in the presence of arrhythmia.

The aim of this work was to propose a novel paradigm that the heart is a communication system where information is shared by electrical wave propagation between cardiomyocytes (figure 1*a*). Under this paradigm, cardiac arrhythmia can be viewed as a state of abnormal production and transmission of information that can be quantified by information theory measures. The specific objective of this work was to develop a framework to quantify cardiac electrical communication during action potential propagation in normal and abnormal heart rhythms. Three major mechanisms of clinically important arrhythmias that could lead to sudden death were considered: *anatomical reentry* (reentry around an anatomically defined circuit [18]), *spiral reentry* (functional reentry without an anatomically defined circuit [19]) and *multiple reentry* (multiple functional reentry circuits in the presence of ongoing wavebreak [20]).

## 2. Material and methods

### 2.1. Models of action potential propagation

To develop the information theory framework, we employed two commonly used mathematical models of cardiac action potential propagation. The first model was a monodomain reaction–diffusion (RD) model that was originally derived by Fitzhugh [21] and Nagumo *et al*. [22] as a simplification of the biophysically based Hodgkin–Huxley equations describing current carrying properties of nerve membranes [23] which was later modified by Rogers & McCulloch [24] to represent cardiac action potential. This model reproduces several physiological properties known to be important in arrhythmogenesis, including slowed conduction velocity (CV) and unidirectional block owing to wavefront curvature [24]. This model was used widely in previous studies [25–31].
2.1
2.2Here, *v* is the excitation variable that can be identified with transmembrane potential, *r* is the recovery variable, *I*_{ex} is the external current [6], and *G _{x}* and

*G*are the conductance in

_{y}*x*- and

*y*-directions on the lattice, respectively. In this study, the lattice was assumed to be isotropic (i.e.

*G*=

_{x}*G*). The model equations were solved using a finite difference method for spatial derivatives and explicit Euler integration for time derivatives assuming Neumann boundary conditions.

_{y}Cellular automata (CA) models have been used to study cardiac action potential propagation in several previous studies [32,33]. The model we used was developed by Alonso Atienza *et al.* [34] who employed realistic restitution properties and the curvature phenomenon. Each cell in the CA can adopt one of the following three states: resting, refractory_{1} and refractory_{2}. Cells in the resting state are relaxed and can be excited, whereas cells in both refractory states are excited. Cells in refractory_{1} can excite neighbouring cells, whereas cells in refractory_{2} cannot. All three states have physiological relevance. The depolarization (or excitation) of a cell *i* is the transition from the resting state into the refractory_{1} state and occurs according to a probabilistic update rule based on two influences [34]: (i) the intrinsic cell excitability (*E*) of *i* that gets higher with the time in the resting state and (ii) the amount of excitation in the neighbouring cells (*Q*):
2.3

where *A _{j}* is the excitation state of cell

*j*(adjacent to

*i*) with a value of 0 for the resting state and a value of 1 for either of the two refractory states, and

*d*is the distance in the lattice between the centre of both cells

_{ij}*i*and

*j*. The transitions from refractory

_{1}to refractory

_{2}(partial repolarization) and from refractory

_{2}to resting (total repolarization) are deterministic. The total time spent in the two refractory states matches the total duration of the action potential (APD). The period in the refractory

_{1}state is equal to 10% of the APD. The parameters

*E*and APD were estimated from restitution curves of the CV and APD, respectively, which depend on the period of the preceding DI [34]. In addition to APD and CV restitution properties, equation (2.3) also reproduces CV slowing in areas with a pronounced wavefront curvature because of the decreased probability of excitation.

### 2.2. Cardiac simulation

For the RD and CA models, the cardiac tissue was simulated as a two-dimensional 128 × 128 isotropic lattice of cells (figure 1*b*). In each cell, the time series of cardiac excitation was computed for 10 s with a sampling frequency of 500 Hz (temporal resolution Δ*t* = 2 ms, figure 1*c*). The sampling frequency and duration of the time series were determined to reflect realistic measurements in human clinical electrophysiology studies [35].

We simulated four different heart rhythms in both models: normal heartbeat, anatomical reentry, spiral reentry and multiple reentry. The latter two are considered to be important mechanisms of cardiac fibrillation, including atrial fibrillation (AF) and ventricular fibrillation (VF) [19,20]. Cardiac simulation was performed using Matlab R2014a (Mathworks, Inc.).

Normal heartbeat was simulated as regular point stimulations (60 beats min^{−1}) originating from the top middle region of the lattice. This pattern of stimulation caused a regular train of curved excitation wavefronts travelling from top to bottom along the vertical axis in both the RD (electronic supplementary material, video S1) and CA models (electronic supplementary material, video S2).

Anatomical reentry is characterized by an electrical wavefront that travels along a preformed anatomical obstacle, most commonly a scar resulting from healed myocardial infarction, and re-excites previously excited tissue. We simulated a cardiac impulse that can rotate around the obstacle, leading to repetitive, rapid excitation of the heart. Anatomical reentry was reproduced in the RD (electronic supplementary material, video S3) and the CA models (electronic supplementary material, video S4) by simulating a non-excitable circular region in the centre of the lattice occupying 20% of the total surface area.

For spiral reentry [19], we simulated a two-dimensional wave of excitation emitted by an organizing source (or *rotor*) of functional reentry whose front is an involute spiral with increasing convex curvature towards the rotation centre [36]. The spiral reentry was generated by a cross-field stimulation protocol [6] in both the RD (electronic supplementary material, video S5) and the CA models (electronic supplementary material, video S6).

Multiple reentry is characterized by multiple independent circuits of functional reentry occurring simultaneously and propagating randomly throughout the cardiac tissue [37]. Wavefronts continuously undergo wavefront–wavetail interactions resulting in wavebreak and generation of new wavefronts [38]. Multiple reentry was reproduced in the RD (electronic supplementary material, video S7) and the CA models (electronic supplementary material, video S8) by a train of random point stimulations in the substrate where the APD is shortened by 40%.

### 2.3. Information measures

For each cell, the time series of cardiac excitation was coarse-grained to 1 when excited (during the APD at 90% repolarization, or APD_{90}) or 0 when resting (figure 1*c*). We treated each cell on the lattice as a time-series process *X*, where at any observation time *t* the process *X* is either excited or resting in that case we define *X _{t}* = 1 or

*X*= 0, respectively.

_{t}Using this framework, we can compute the Shannon entropy *H* of each time-series process *X*where *p*(*x*) denotes the probability density function of the time series generated by *X*. This quantifies the average uncertainty of whether a single cell is excited (*x* = 1) or resting (*x* = 0) over each cell's time history [39].

Mutual information *I*(*X*; *Y*) is a measure of the reduction in uncertainty of the time-series process *X* owing to the information gained from knowing the time-series process *Y*; hence, this quantity is commonly viewed as the information shared between *X* and *Y* [39]. Therefore, by computing the mutual information, we can receive insights into which cells share information and how much. Formally, the mutual information between the time-series processes, in this case cells, *X* and *Y* is
2.4
2.5where *p*(*x*, *y*) and *H*(*X*, *Y*) denote the joint probability density function and the joint entropy of *X* and *Y*, respectively (figure 1*a*).
2.6
2.7
To understand the spatial profiles of information sharing between cardiomyocytes, mutual information was computed between each of five representative cells (green circles in figure 1*d* and all other cells in the two-dimensional lattice taken individually). These representative cells were defined to be in the left-upper quadrant (`32,32`), the right-upper quadrant (`32,96`), the centre (`64,64`), the left-lower quadrant (`96,32`) and the right-lower quadrant (`96,96`). These points were chosen to avoid artefacts generated by the boundary conditions and point stimulation in the RD model as discussed in the Results section. Custom programs in Python were used to compute information measures.

## 3. Results

### 3.1. Normal heartbeat

In the RD model, electrical wavefronts regularly swept the lattice from top to bottom (electronic supplementary material, video S1 and figure 2*a*, top row). The entropy was relatively lower at the lattice borders and higher at the site of stimulation, but these were artefacts of the boundary conditions and point stimulation, respectively (figure 2*b*, top row). Otherwise, entropy was homogeneous across the lattice [0.68 (mean) ± 0.03 (s.d.) bits]. Mutual information showed a spatially heterogeneous information sharing between cells (figure 2*c*, top row) despite the assumed isotropic structure and homogeneous electrical properties of the lattice. For example, figure 2*c*3 (top row) shows that the cell in the centre of the lattice shares a high amount of information with the cells on the same electrical wavefront (yellow band) generated by the heartbeats. It also shows little information sharing with the cells preceding and following the wavefront (light blue bands surrounding the yellow band). This is clearly shown in the profile of mutual information (figure 2*d*, top row) along the vertical broken line in figure 2*c*3 (top row). Mutual information (red line) reached its peak (0.69 bits) at the centre (profile position 64), where mutual information is equal to entropy, because the mutual information of an entity with itself is equal to entropy (equation (2.5)). Mutual information fell off sharply from the centre and reached the minimum (0 bits) on both sides (profile position 41 and 86) before slightly rising to approximately 0.06 bits on both ends. These findings indicate that the cells share a high amount of information when they are in phase with the cardiac excitation and share little information when they are out of phase. We formed an analytical framework which corroborates these numerical results (electronic supplementary material, appendix S1).

In the CA model (electronic supplementary material, video S2 and figure 2*a*, bottom row), the electrical wavefront was more irregular and unstable than that of the RD model. Entropy was homogeneous across the lattice, but was lower (0.52 ± 0.00 bits) than that of the RD model (figure 2*b*, bottom row). This difference results from the fact that the CA model had a longer resting state (electronic supplementary material, video S2), making it more biased towards the resting state than the RD model (electronic supplementary material, video S1). There was no qualitative difference in information sharing between the RD model (figure 2, top row) and the CA model (figure 2, bottom row). However, information sharing was lower in the CA model, because the CA model is inherently probabilistic and less reproducible than the RD model. The profile of mutual information (figure 2*d*) along the vertical broken line in figure 2*c*3 also shows a qualitatively similar but lower mutual information in the CA model (figure 2*d*) relative to the RD model (figure 2*d*).

### 3.2. Anatomical reentry

In both the RD (electronic supplementary material, video S3) and the CA models (electronic supplementary material, video S4), the entropy of the cells within the circular non-excitable region was zero, because these cells were always in the resting state (figure 3*a,b*). The entropy was roughly homogeneous in other regions of the lattice. The average entropy was 0.32 ± 0.17 bits in the RD model and 0.69 ± 0.33 bits in the CA model (figure 3*b*). The difference resulted from the longer wavelength in the CA model owing to the more convex curvature than the RD model.

Overall, both models showed a similar spatial pattern of mutual information (figure 3*c*). Information sharing between cells was spatially heterogeneous, but showed rotational symmetry about the non-excitable region in the centre. For example, figure 3*c*1 shows that the cell in the left upper quadrant of the lattice shares a high amount of information with the cells on the same electrical wavefront (the orange band in the RD model and the yellow band in the CA model). Information sharing in the cells on the three other quadrants (figures 3*c*2, 4 and 5) was rotationally symmetric with that of figure 3*c*1. Importantly, there is no information sharing between the cell within the circular non-excitable region and any other cells in the lattice (figure 3*c*3). This is logical from both the standpoints of electrophysiology and information theory. Of note, similar to normal heartbeat (figure 2*c*), both models also showed little information sharing with the cells that preceded and followed the region of a high amount of shared information (light blue bands before and after the yellow band). These findings indicate that, similar to normal heartbeat, the cells share a high amount of information when they are in phase with cardiac excitation, and share little information when they are out of phase.

### 3.3. Spiral reentry

In the RD model (electronic supplementary material, video S5 and figure 4*a*, top row), a single rotor with a resultant spiral reentry was simulated in the lattice. The entropy of each individual cell in this simulation shows an important finding with potential clinical significance. The left lower quadrant and both upper quadrants of the lattice exhibited homogeneous entropy (0.74 ± 0.04 bits), except for the borders of the lattice which is an artefact of boundary conditions (figure 4*b*, top row). The red region in the right lower quadrant represents the higher entropy near the spiral tip (rotor) caused by the CV slowing near the rotor owing to a pronounced wavefront curvature (figure 4*b*, top row). This slowing of CV effectively caused longer cardiac excitations, making the cells in this region more biased towards the excited state than the rest of the cells in the lattice not directly affected by the rotor, boosting the entropy in this region. Within the red region is an L-shaped, light green, bead-like structure representing the lower entropy in the path of the drifting core of the spiral reentry around which the rotor revolved. This result suggests that the entropy of individual cells may be used to aid in localizing the core of spiral reentry for therapeutic purposes. Of note, the L-shaped path of the drifting core in this model was artificially determined by the boundary condition of the model.

While the entropy of individual cells provided very important findings, the average entropy over all cells was less informative. In fact, the average entropy between normal heartbeat and spiral reentry in the RD model was very similar (0.68 ± 0.03 versus 0.74 ± 0.04 bits, respectively). This suggests that the spatial profiles of entropy are more useful in highlighting the difference in dynamics than the aggregate information over the entire tissue, as averaging effectively filters out the important features of the entropy landscape.

Mutual information was far more spatially heterogeneous than anticipated from the electrical wave propagation. For example, figure 4*c*3 (top row) shows that the high level of information sharing in the central region of the lattice quickly faded as distance between cells increased. This reflects the fact that the cells lying along the same electrical wavefront changed over time owing to the drifting core. This resulted in the smaller region of high information sharing in the spiral reentry than in the normal heartbeat. This finding indicates that mutual information can sensitively detect regional heterogeneity of cardiac excitation in spiral reentry, which is not apparent from electrical wave propagation. Of note, information sharing in the right lower quadrant was limited to a focal region without a spiral tail (figure 4*c*5, top row). This is because the cell in the right lower quadrant (green circle) happened to lie on the path of the drifting core, which coincided with a void of cardiac excitation.

In the CA model (electronic supplementary material, video S6 and figure 4*a*, bottom row), the electrical wavefront was inherently more irregular than that of the RD model. The entropy was roughly homogeneous across the lattice and close to 1 (0.97 ± 0.00 bits; figure 4*b*, bottom row). The higher entropy of the CA model compared with the RD model was, as in the anatomical reentry, caused by the longer wavelength in the CA model owing to the more convex curvature than the RD model (figure 4*a*, bottom row). Unlike the RD model, the entropy in the CA model did not show the path of the drifting core, indicating that the drift trajectory was much more random compared with that of the RD model with respect to the time frame (10 s) of data acquisition. Information sharing between cells was spatially heterogeneous (figure 4*c*,*d*) because of the regional heterogeneity of cardiac excitation in spiral reentry. Overall, the CA model showed lower information sharing than the RD model owing to the probabilistic nature of the model.

### 3.4. Multiple reentry

In both the RD (electronic supplementary material, video S7) and the CA models (electronic supplementary material, video S8), the entropy was homogeneous across the cell lattice (figure 5*a*). The average entropy was 0.88 ± 0.03 bits in the RD model and 0.97 ± 0.00 bits in the CA model (figure 5*b*). This indicates that excited and resting states are almost equally distributed throughout the time series of all the cells, yielding a high uncertainty and homogeneous entropy. The spatial profiles of information sharing for both the RD and the CA models were similar to that of spiral reentry (figure 5*c*), except the fact that the underlying structure was much less organized owing to the random nature of multiple reentry. Information sharing was low except for a small region in the immediate neighbourhood of the cell in which mutual information was measured. Outside this small region information sharing steeply fell off to near zero (figure 5*d*). The near-zero mutual information indicates that the cells almost completely lost synchrony during multiple reentry; that is, individual cardiomyocytes got excited independently from each other and did not share information with cells beyond their immediate neighbourhood.

## 4. Discussion

### 4.1. Summary of the findings

By treating the heart as an electrical communication system, we demonstrated quantitatively that information sharing between cardiomyocytes on an isotropic lattice structure is spatially heterogeneous. This finding was unexpected from the traditional concept of the heart as a functional syncytium sharing electrical information via gap junctions, where one might mistakenly assume that information sharing would be homogeneous along the electrical wavefront. We also found that entropy can be significantly different between heart rhythms with electrically similar spatial patterns (figures 2*b*, 4*b* and 3*b*). These findings indicate that metrics from information theory can quantitatively assess the communication processes within the heart which are not obvious from conventional electrocardiographic metrics such as sequences of electrical excitation. In addition, our results show that cardiac arrhythmia significantly impacts electrical communication within the heart.

### 4.2. Mutual information to quantify communication within the heart

Analysis of dynamical multivariate datasets over the dimensions of time and physical space is commonly encountered in the investigation of complex systems [40]. The study of cardiac arrhythmia, particularly cardiac fibrillation, is no exception. For example, a number of measures to quantify spatial complexity of VF have been proposed, including the correlation length [41], the multiplicity index [42] and Karhunen–Loève decomposition [43]. The main focus of interest in these studies was to quantify the determinism and the predictability of the time series over physical space.

Our study is different from these previous studies in two aspects. First, our focus of interest was to quantify communication within the heart. Of central importance to the understanding of complex systems is connectivity, or the presence of dynamical interactions between spatially distinct locations within the system. Knowledge about connectivity in a system, whether anatomical or functional, further facilitates the fundamental understanding of the system because it addresses an important aspect of the functional interdependency between each component of the system. Our results indicate that information theory metrics can quantitatively assess electrical communication processes among cardiomyocytes during normal heartbeat and complex arrhythmias beyond electrocardiographic measures, conferring validity to the paradigm of the heart as a communication system. Second, we used mutual information to perform spatial profiling of different cardiac arrhythmias. Correlations within multivariate time series can be described by measures such as linear and nonlinear correlation functions. However, mutual information has attracted considerable attention recently because it promises a very general quantification of statistical dependence [44]. In addition, previous studies measured the spatial profile of mutual information during VF [45], which could include both spiral reentry and multiple reentry, because it is challenging to distinguish one from the other in experimental settings. Therefore, the spatial profile of spiral reentry and multiple reentry was not clearly delineated. Our result showed that the spatial profile and the underlying structure of mutual information during each arrhythmia are clearly different (figures 4 and 5). This suggests that information theory metrics could help distinguish one rhythm from another by quantifying communication within the heart.

The underlying mechanism of perturbation of information transfer during arrhythmia remains unclear. Information sharing during VF seems to be directly affected by the anisotropy of myofibre orientation and cell-to-cell coupling [45]. However, the spatial profile of membrane potential during VF has no consistent relationship with that of intracellular calcium dynamics [46]. This may suggest a contribution of non-voltage-gated intracellular calcium release in perturbation of information transfer by increasing the local complex interactions between calcium dynamics and membrane potential. Clearly, further studies will be needed to investigate the mechanistic basis of the paradigm of the heart as a communication system.

### 4.3. Clinical implications

Recently, targeted elimination of the rotor (phase singularity) of spiral reentry has been shown to result in sustained termination of AF [47]. As a result, spatial localization of the rotor has attracted substantial attention in clinical cardiac electrophysiology. Entropy has been quantified to identify the location of the rotor of spiral reentry from the bipolar electrograms by creating the probability density function based on the amplitude of the signal [48]. However, the accuracy of this metric was not clear, because it has consistent correlation with complex fractionated electrograms [49], which was found to bear no spatial relevance to spiral reentry [50]. Therefore, the knowledge of the spatial profile of entropy for spiral reentry was lacking.

Our result in the RD model clearly showed that the region of rotor drift has high entropy (red region, right lower quadrant in figure 4*b*, top row), which is consistent with previous studies [48]. Moreover, what was most striking was the fact that entropy can localize the path of the drifting core of spiral reentry (L-shaped, light green, bead-like structure, right lower quadrant in figure 4*b*, top row) because of the low entropy of the spiral core. This makes electrophysiological sense, because the cardiomyocytes within the spiral core are almost constantly depolarized [51], making the probability density-biased towards 1. Therefore, our result showed a critically important fact that entropy can spatially localize the core (low entropy) within a larger region of the drifting rotor (high entropy). However, because a similar structure of the spatial profile could not be identified in the CA model (figure 4*b*, bottom row), localization of the drifting core may require a spatially stable spiral reentry with adequately slow drift. Although these preliminary findings need to be confirmed in experimental models of cardiac fibrillation, they illustrate the potential clinical utility of information theory applied to cardiac electrophysiology.

### 4.4. Limitations

There are two limitations that should be considered before our results can be translated to human patients. First, the cardiac tissue was assumed to be a two-dimensional, isotropic and homogeneous lattice, whereas real cardiac tissue is three-dimensional, anisotropic and heterogeneous owing to the intricately woven myofibre structure [52] and regional heterogeneity [53]. These tissue properties may contribute critically to the generation of cardiac arrhythmia [54]. However, the main focus of this work was to prove the concept that quantitative analysis of electrical communication during existing cardiac arrhythmia could yield clinically relevant results. We used two widely accepted models of action potential propagation in cardiac tissue to reproduce a variety of heart rhythms that captured important features of clinically representative arrhythmias. Therefore, we believe that these model assumptions were acceptable within the scope of this work. Second, our computation of information theory metrics did not incorporate conduction delay of electrical current to travel from one cell to another, because the time series of the entire lattice was acquired simultaneously. This is because the conduction delay within the small two-dimensional lattice would be negligibly small relative to the acquisition period of 10 s. However, this assumption may have underestimated the true amount of information sharing between heart cells, because the standard definition of mutual information does not include shared information that is delayed in time. This leaves open the potential for even more clinically useful results by considering generalizations of the metric that explicitly account for this conduction delay.

## 5. Conclusion

Information theory metrics can quantitatively assess electrical communication processes among cardiomyocytes during normal heartbeat and complex arrhythmias beyond electrocardiographic measures. Further, entropy may have a clinical application in the localization and elimination of spiral reentry cores. These results suggest that the heart as a communication system is more complex than the traditional concept of functional syncytium sharing electrical information via gap junctions. We believe that this new paradigm provides a new set of tools for the systems-approach to the heart as a complex system [55].

## Acknowledgements

This work was conducted as a research project during the Complex Systems Summer School at the Santa Fe Institute (SFI), NM, USA. The authors thank all the staff at SFI, particularly Sander Bais, Juniper Lovato, and John-Paul Gonzales. The authors also thank Felipe Alonso Atienza, Ferney Beltrán-Molina, Jesús Requena Carrión and Peter Hammer for generously providing the source code for the models described in the paper. The authors also thank Simon DeDeo and Nix Barnett for valuable input.

- Received October 31, 2014.
- Accepted February 11, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.