Classifying dynamic transitions in high dimensional neural mass models: A random forest approach

<div><p>Neural mass models (NMMs) are increasingly used to uncover the large-scale mechanisms of brain rhythms in health and disease. The dynamics of these models is dependent upon the choice of parameters, and therefore it is crucial to be able to understand how dynamics change when parameters are varied. Despite being considered low dimensional in comparison to micro-scale, neuronal network models, with regards to understanding the relationship between parameters and dynamics, NMMs are still prohibitively high dimensional for classical approaches such as numerical continuation. Therefore, we need alternative methods to characterise dynamics of NMMs in high dimensional parameter spaces. Here, we introduce a statistical framework that enables the efficient exploration of the relationship between model parameters and selected features of the simulated, emergent model dynamics of NMMs. We combine the classical machine learning approaches of trees and random forests to enable studying the effect that varying multiple parameters has on the dynamics of a model. The method proceeds by using simulations to transform the mathematical model into a database. This database is then used to partition parameter space with respect to dynamic features of interest, using random forests. This allows us to rapidly explore dynamics in high dimensional parameter space, capture the approximate location of qualitative transitions in dynamics and assess the relative importance of all parameters in the model in all dimensions simultaneously. We apply this method to a commonly used NMM in the context of transitions to seizure dynamics. We find that the inhibitory sub-system is most crucial for the generation of seizure dynamics, confirm and expand previous findings regarding the ratio of excitation and inhibition, and demonstrate that previously overlooked parameters can have a significant impact on model dynamics. We advocate the use of this method in future to constrain high dimensional parameter spaces enabling more efficient, person-specific, model calibration.</p></div

Common dynamic patterns observed in the Wendling model.

<p>Non-steady-state solutions are split into two categories: oscillations and (poly)spike-wave dynamics. Oscillations are cycles with one peak per period delineated by frequency into the five classical clinical bands: gamma (30-60Hz); beta (13-30Hz); alpha (8-12Hz); theta (4-8Hz); and delta (0-4Hz). (Poly)spike-wave dynamics are cycles with one or more spikes per period, riding on an oscillation of between 2 and 8 Hz with a mean of 4 Hz. For the sake of clarity amplitudes do not have uniform y-axis scales.</p

Bivariate joint distribution of the likelihood of steady state (lower triangle) or seizure dynamics (upper triangle).

<p>Each subfigure is a projection of the parameter space over two parameters, and the colour indicates the likelihood of finding a particular type of dynamics (seizure or steady state) as per the colourbar. For example, the subfigure in the second column on the first row (encircled and labelled (i)) maps the likelihood of finding seizure dynamics over different values of <i>B</i>(x-axis) and <i>A</i>(y-axis), given variations in all other parameters. In the upper triangle, yellow indicates high likelihood of observing seizure dynamics, whereas blue indicates low likelihood of observing seizure dynamics. In the lower triangle, red indicates high likelihood of observing steady state dynamics, whereas blue indicates low likelihood of observing steady state dynamics. Each subfigure was computed using <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1006009#pcbi.1006009.e015" target="_blank">eq 12</a> with 20 × 20 bins over the parameter ranges provided in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1006009#pcbi.1006009.t002" target="_blank">Table 2</a>. Upper triangle: Specific combinations of parameters can lead to manifolds with a high likelihood of seizure dynamics (see for example the linear relationships between <i>A</i> and <i>B</i> in the encircled subfigure (i) and <i>a</i> and <i>b</i> in the encircled subfigure (ii)). Lower triangle: one can observe that small values of the parameters <i>A</i> or <i>C</i> guarantee a steady state (see for example the encircled subfigures (iii) and (iv)).</p

A tree representing how the extended parameter space (incorporating two additional ‘ratio parameters’ <i>r</i>_{<i>a</i>/<i>b</i>} and <i>r</i>_{<i>A</i>/<i>B</i>}) is split dependent on the presence or absence of seizure dynamics.

<p>The root region is at the top of the figure and represents the total parameter space, while the leaves are at the bottom. The upper label of each region indicates the size of the parameter space represented in this region. The lower label indicates the percentage of all parameter combinations that result in seizure dynamics. One can see that the ratios have an important role in splitting the parameter space. Values are given to two digits precision.</p

The importance of parameters as determined by normalised variable importance (<i>NVI</i>) averaged over a random forest of 100 trees.

<p>The ratios <i>r</i>_{<i>A</i>/<i>B</i>} and <i>r</i>_{<i>a</i>/<i>b</i>} have been added as additional parameters. Four characteristics of interest are considered: the switch between steady state and non-steady state, the amplitude of cycles, the frequency of cycles and the switch between any activity (mainly steady state) and seizure dynamics. A value of 1 signifies the parameter with greatest importance for observing the feature of interest. A value of 0 implies a parameter has no control over observing a feature of interest.</p

A tree representing how parameter space is split dependent on the presence or absence of seizure dynamics.

<p>The root region is at the top of the figure and represents 100% of the parameter space while the leaves are at the bottom. The upper label of each regions indicates the size of the parameter space represented in this region. The lower label indicates the percentage of all parameter combinations that result in seizure dynamics. The colour indicates the density of seizure dynamics in the given region. The parameters <i>A</i>, <i>B</i>, <i>a</i> and <i>b</i> are the most important parameters because they efficiently split the parameter space into subspaces with high or low likelihood of seizure dynamics. Some parameters such as <i>P</i> and <i>C</i> do not appear in this small tree, however they can appear in a more complex tree (see supplementary material, <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1006009#pcbi.1006009.s002" target="_blank">S1 Fig</a>). Values are given to two digits precision.</p

The range of considered parameter space of the Wendling model.

<p>Details of the reference used to define the minimum and maximum value of each parameter is included. Chosen ranges were constrained either by experiments (e.g. <i>a</i> and <i>b</i>) or the widest range described in theoretical studies (e.g. <i>P</i> and <i>C</i>).</p