# Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers

David R Penas, Meysam Hashemi, Viktor K Jirsa, Julio R Banga

► To cite this version:

David R Penas, Meysam Hashemi, Viktor K Jirsa, Julio R Banga. Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers. PLoS Computational Biology, 2024, 20, 10.1371/journal.pcbi.1011642. [https://doi.org/10.1371/journal.pcbi.1011642](https://doi.org/10.1371/journal.pcbi.1011642)

HAL Id: hal-04890371

https://hal.science/hal-04890371v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers

David R. Penas¹*, Meysam Hashemi²*, Viktor K. Jirsa²*, Julio R. Banga¹*

¹ Computational Biology Lab, MBG-CSIC (Spanish National Research Council), Pontevedra, Spain, ² Aix Marseille Univ, INSERM, INS, Inst Neurosci Syst, Marseille, France

* These authors contributed equally to this work.

* viktor.jirsa@univ-amu.fr (VKJ); j.r.banga@csic.es (JRB)

<figure><img src="image_2.png" /></figure>

## Abstract

The Virtual Epileptic Patient (VEP) refers to a computer-based representation of a patient with epilepsy that combines personalized anatomical data with dynamical models of abnormal brain activities. It is capable of generating spatio-temporal seizure patterns that resemble those recorded with invasive methods such as stereoelectro EEG data, allowing for the evaluation of clinical hypotheses before planning surgery. This study highlights the effectiveness of calibrating VEP models using a global optimization approach. The approach utilizes SaCeSS, a cooperative metaheuristic algorithm capable of parallel computation, to yield high-quality solutions without requiring excessive computational time. Through extensive benchmarking on synthetic data, our proposal successfully solved a set of different configurations of VEP models, demonstrating better scalability and superior performance against other parallel solvers. These results were further enhanced using a Bayesian optimization framework for hyperparameter tuning, with significant gains in terms of both accuracy and computational cost. Additionally, we added a scalable uncertainty quantification phase after model calibration, and used it to assess the variability in estimated parameters across different problems. Overall, this study has the potential to improve the estimation of pathological brain areas in drug-resistant epilepsy, thereby to inform the clinical decision-making process.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Penas DR, Hashemi M, Jirsa VK, Banga JR (2024) Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers. PLoS Comput Biol 20(7): e1011642. https://doi.org/10.1371/journal.pcbi.1011642

**Editor:** Maxim Bazhenov, University of California San Diego, UNITED STATES OF AMERICA

**Received:** October 31, 2023

**Accepted:** June 18, 2024

**Published:** July 11, 2024

**Copyright:** © 2024 Penas et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The patient data sets cannot be made publicly available due to data protection concerns regarding potentially identifying and sensitive patient information. Interested researchers may access the data sets by contacting Clinical Data Manager Aurelie Ponz (aurelie.ponz@univ-amu.fr) at the Institut de Neurosciences des Systemes, Aix-Marseille Universite. Code availability: Detailed results and the main source codes needed to reproduce them are available on Zenodo https://doi.org/10.5281/zenodo.10057788. For the most recent version of

## Author summary

Motivated by the problem of parameter estimation in a set of whole-brain network models of epilepsy (of increasing complexity), this study addresses the question of choosing a robust global optimization solver that can be accelerated by exploiting parallelism in different infrastructures, from desktop workstations to supercomputers. By leveraging data-driven techniques with robust cooperative global optimization methods, we aim to achieve accurate parameter estimation with reduced reliance on prior information. This is due to the dependency of Bayesian inference on the level of information in the prior, while this approach allows us to quantify uncertainty in the absence of any prior knowledge effectively. In this work, we construct an efficient and accurate method to perform
the code, visit GitHub at https://github.com/ins-
amu/BVEP/tree/master/Optimization.

**Funding:** This work was funded by the Spanish Agencia Estatal de Investigación MCIN/AEI/ 10.13039/501100011033 grant PID2020- 117271RB-C22 BIODYNAMICS (https://www.aei. gob.es/) to JRB, by the French National Research Agency (ANR) as part of the second "Investissements d'Avenir" program, ANR-17- RHUS-0004, EPINOV (https://anr.fr) and ANR-22- PESN-0012 (France 2030 program) to VJ and MH, by the European Union's Horizon Europe Programme under the Specific Grant Agreement No. 101147319 EBRAINS 2.0 Project (https://ec. europa.eu/programmes/horizon2020) to VJ. The funders did not play any role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

parameter estimation and uncertainty quantification for the whole-brain models, and we use it to infer the brain regional epileptogenicity from source and sensor level data. Of specific interest is the ability of our method to produce inference for high-dimensional state-space models governed by deterministic, stochastic, well-behaved, and stiff differential equations, using only partial observations and sparse encoding from system states to the observation.

Introduction

Achieving maximum accuracy in predictive modeling poses a significant challenge in inferring critical information from recording data, such as characterizing pathological brain networks from (stereoelectro) EEG data. Optimizing parameters efficiently is therefore essential for gaining insights into underlying mechanisms of the data-generating process in pathological conditions, such as the dynamics and connectivity of abnormal brain networks.

The quantitative approaches have accompanied a major conceptual and theoretical evolu-
tion in epilepsy surgery, moving beyond the concept of an epileptogenic focus and proposing
several definitions of different zones (e.g., the epileptogenic zone, the seizure onset zone, and
the irritative zone [1-5]). In partial epilepsy, seizures originate from a network of hyperexcit-
able regions referred to as the epileptogenic zone (EZ; [6]) and then propagate to a secondary
connected network, the so-called propagation zone (PZ; [7]). However, defining the EZ based
on stereoelectro EEG (SEEG) can be difficult, as the patterns of seizure onset and early spread
may exhibit complex organization, influenced by network dynamics [4, 8]. For instance, it has
been reported that the seizure onset zone (SOZ) is not always a good marker of the EZ [9].
Given that epileptogenic regions often originate beyond a single epileptic focus, a network of
brain regions that is necessary and sufficient for initiation of seizures, the removal of which
guarantees the complete abolition of seizures [5, 10]. In this study focusing on model-based
inference using in-silico data, we assume that the EZ is equivalent to the SOZ. Nevertheless,
the success of surgical interventions for drug-resistant patients critically depends on the preci-
sion and reliability of the initial hypotheses, e.g., the spatial map of EZ/PZ as an identification
of the seizure organization [5, 10-15].

Virtual Epileptic Patient (VEP; [10, 16]) is a digital modeling approach that integrates mathematical modeling of abnormal neural activity with patient-specific anatomical data to predict the brain network involved in seizure generation and propagation in individuals. By accurately capturing diverse seizure dynamics and generating computer simulations resem- bling intracranial EEG recordings, this technique offers a versatile platform to optimize the surgical strategy and to aid in clinical decision-making [5, 10, 17, 18]. The VEP is a model- based approach that relies on estimating the parameters in a high-dimensional state-space representation to accurately identify the network of EZ/PZ in the brain. Besides the issues of sparsity, stochasticity, and scalability, the reliability of prediction on the EZ/PZ is challenging due to the nontrivial effects of brain networks, the non-linearity involved in the spatiotemporal organization of the brain, and the uncertainty in model components.

Several studies [5, 10, 13, 17, 19–22] have demonstrated that advancements in more accu-
rate estimation of the VEP parameters have the potential to more informed clinical decision-
making and optimize surgical strategies. This motivation drove us to benchmark various paral-
lel global search optimization techniques that eliminate the need for strong initial assumptions,
such as informative prior information in a Bayesian setup. Although clinical information or
knowledge extracted from multimodal data integrated into the prior can significantly enhance
the accuracy of EZ/PZ estimation [13, 20], the use of global optimization methods becomes
crucial in scenarios where such prior knowledge is lacking. Global search optimization meth-
ods offer several advantages, such as an scalability and flexibility in objective functions and
constraints, enhanced exploration-exploitation trade-off, robustness to initial conditions, and
the capability to deal with non-convex and multimodal problems. By leveraging the advantages
of global search optimization methods with high-performance computing (HPC) infrastruc-
tures, we can achieve faster and more accurate parameter estimation, thereby enhancing the
overall efficacy of our approach.

**Materials and methods**

**Individual patient data**

For this study, we use secondary data from a selected patient (a 35 year-old male) initially diag-
nosed with left temporal lobe epilepsy (Histopathology: Gliosis, Surgical procedure: resection,
Surgical outcome: seizure free, Engel score I). The patient underwent comprehensive presurgi-
cal evaluation, including clinical history, neurological examination, neuropsychological test-
ing, structural and diffusion MRI scanning, Stereotactic-EEG (SEEG) recordings along with
video monitoring as previously described in [7, 13]. The evaluation included non-invasive
T1-weighted imaging (MPRAGE sequence, repetition time = 1900 ms, echo time = 2.19 ms,
1.0 x 1.0 x 1.0 mm, 208 slices) and diffusion MRI images (DTI-MR sequence, angular gradient
set of 64 directions, repetition time = 10.7 s, echo time = 95 ms, 2.0 x 2.0 x 2.0 mm, 70 slices, b-
weighting of 1000 smm⁻²). The images were acquired on a Siemens Magnetom Verio 3T MR-
scanner.

**Network anatomy**

The structural connectome was built using the TVB-specific reconstruction pipeline with gen-
erally available neuroimaging software, as described in [17, 23]. First, the command *recon-all*
from the Freesurfer package [24] in version v6.0.0 was used to reconstruct and parcellate the
brain anatomy from T1-weighted images. Then, the T1-weighted images were coregistered
with the diffusion-weighted images using the linear registration tool *flirt* [25] from the FSL
package in version 6.0, with the correlation ratio cost function and 12 degrees of freedom.

The MRtrix package version 0.3.15 was then used for tractography. The fiber orientation
distributions were estimated from diffusion-weighted images using spherical deconvolution
[26] by the *dwi2fod* tool, with the response function estimated by the *dwi2response* tool using
the *tournier* algorithm [27]. Next, we used the *tckgen* tool, employing the probabilistic tracto-
graphy algorithm iFOD2 [28] to generate 15 million fiber tracts. Finally, with the generated
fiber tracts and the regions defined by the brain parcellation, the connectome matrix is built
by counting the fibers connecting all regions. Using the *tck2connectome* tool and the Desikan-
Killiany parcellation [29], the patient's brain is divided into 68 cortical regions and 16 subcorti-
cal structures. See S1 File for the label names and indices of the sub-divided brain regions. The
connectome was normalized so that the maximum value is equal to one (see S1 File).

**Virtual Epileptic Patient (VEP) model.** In the process of building a personalized brain network model, the brain regions (network nodes) are defined using a parcellation scheme, and a set of dynamical equations, known as the neural mass model, is placed at each network node to generate the regional brain activity [16, 30]. Neural masses are commonly used to model the collective behavior of populations of neurons in the brain, such as firing rates, capturing macroscopic dynamics and interactions rather than focusing on individual neuron behavior [31–34]. They have demonstrated efficiency in capturing the main features of brain functional behaviors in a single computational framework by accounting for interactions
among brain regions in healthy and diverse pathological conditions [10, 21, 30, 35–38]. Taking a data-driven approach that integrates subject-specific brain anatomy, the network edges are subsequently encoded using a personalized structural connectivity matrix derived from non-invasive imaging data, such as diffusion magnetic resonance imaging (dMRI), for an individual subject [16, 17]. Moreover, the anatomical data imposes a constraint on simulated data, wherein the salient aspects of the dynamics emerged at each brain region, allowing the hidden state dynamics to be inferred from the data. In the Virtual Epileptic Patient (VEP model; [16]), a personalized brain network model of epilepsy spread, the dynamics of brain network nodes are governed by the so-called Epileptor model [39]. The Epileptor is a general description of epileptic seizures (in humans, mice, rats, and zebrafishes), which contains the complete taxonomy of system bifurcations to realistically reproduce the dynamics of onset, progression, and offset of seizure-like events [40, 41]. The full Epileptor model consists of five state variables that couple two oscillatory dynamical systems operating on three distinct time scales. At the fastest time scale, an oscillatory dynamical system accounts for fast discharges during ictal seizure states, whereas on the intermediate time scale, another system represents the slow spike-and-wave oscillations. On the slowest time scale, a permittivity state variable is responsible for the transition between interictal and ictal states. The permittivity variable represents the slow-evolving extracellular processes that occur during epileptiform activity, such as levels of ions, oxygen, and energy metabolism. Depending on its values, the dynamics of Epileptor may drive it into or out of a seizure, which accounts for its bi-stable behavior (for more details, see [39]).

Motivated by Synergetic theory [42, 43] and assuming a time-scale separation ($\tau \gg 1$), the fast variables swiftly converge onto the slow manifold, which is governed by the dynamics of the slow variable. This adiabatic approximation [44, 45] leads to the 2D reduction of VEP model, given by:

$$ 
\begin{aligned}
\dot{x}_i &= 1 - x_i^3 - 2x_i^2 - z_i + I_i \\
\dot{z}_i &= \frac{1}{\tau} (4(x_i - \eta_i) - z_i - K \sum_{j=1}^{N} C_{ij}(x_j - x_i)),
\end{aligned}
\quad (1) $$

where $x_i$ and $z_i$ indicate the fast and slow variables corresponding to i-th brain region, respectively. The parameter $I = 3.1$ represents the flow of electric current, and $\tau$ scales the length of the seizure. The degree of epileptogenicity at each brain region is represented by the value of the excitability parameter $\eta_i$ (Hopf bifurcation parameter). The network nodes are coupled by a linear diffuse approximation of permittivity coupling through $K \sum_{j=1}^{N} C_{ij}(x_j - x_i)$, which includes a global scaling factor $K$ on the subject's connectome $C_{ij}$. In this paper, this set of VEP equations is solved using a Euler integration scheme with a step size of 0.1 msec. The initial conditions for the $x$ and $z$ variables at each brain region were selected as (-2.5, 3.5), respectively.

**Spatial Map of Epileptogenicity.** In the VEP model, the occurrence of seizures in each brain region depends on its individual excitability (node dynamics) and its connectivity to other regions (network coupling). The excitability parameter, denoted by $\eta$, regulates the level of tissue excitability, and its spatial distribution is the focus of parameter estimation. The brain regions can be classified into three main types based on their excitability parameter:

* **Epileptogenic Zone (EZ):** If $\eta > \eta_c$, the brain region will autonomously trigger seizures, which are responsible for initiating and organizing the early stages of epileptic activity. The Epileptor model exhibits an unstable fixed point in these regions, allowing seizures to occur independently of network effects.
* Propagation Zone (PZ): If $\eta_c - \Delta\eta < \eta < \eta_c$, the brain region does not autonomously trigger seizures. However, it can be recruited in seizure propagation since its equilibrium state is close to the critical value. In these regions, a supercritical Andronov-Hopf bifurcation occurs at $\eta = \eta_c$, triggering seizure onset when a sufficiently large external input is present. Otherwise, the Epileptor model remains in a stable equilibrium state.

* Healthy Zone (HZ): If $\eta < \eta_c - \Delta\eta$, the brain region remains seizure-free. In these regions, all trajectories in the phase-plane of the Epileptor model converge to a single stable fixed point, indicating a healthy (non-epileptic) state.

Based on the above dynamical properties, the spatial distribution of epileptogenicity across different brain regions is determined by estimated heterogeneity in excitability parameters: EZs exhibiting high excitability, while PZs have lower excitability values, and very low values of excitability characterize all other regions as HZs.

Note that having an intermediate excitability value (i.e., close to the bifurcation value) does not guarantee recruitment into the seizure propagation. Seizure recruitment is governed by various factors, including structural connectivity, network coupling, and brain state dependence on noise, which all play a crucial role in determining the extent of propagation within the brain network [13, 17]. The linear stability analysis indicates that, in the absence of coupling ($K=0$), the isolated nodes exhibit a bifurcation at the critical value $\eta_c = -2.05$ [46], and we set $\Delta\eta = 1.0$ [13, 17].

## Simulated Stereotactic-EEG (SEEG) data

Simulated Stereotactic-EEG (SEEG) implantation produces data to be used in the building and validation of VEP models. This invasive method is used in clinical situations for patients with drug-resistant epilepsy to determine the focal location of epileptic seizures [2, 5, 6, 10]. The implanted SEEG electrodes record the local field potential generated by the neuronal tissue in its vicinity. The gain matrix (also known as the lead-field or projection matrix) maps the source activity to the measurable sensor signals. Each sensor collects the source signals in its proximity, weighted by the distance and orientation of sources. To model the SEEG signals, here we assume an exponential relation between the source activities and the measurable signals at the sensors:

$$S_i(t) = \sum_{j=1}^{N_n} G_{ij} \exp(x_j(t)), \quad (2)$$

where $S_i(t)$ is the SEEG signal at sensor $i \in \{1, 2, \dots, N_s\}$ with $N_s$ the number of channels (sensors), $x_j(t)$ is the source activity (given by fast variable in Eq 1) in region $j \in \{1, 2, \dots, N_n\}$ with $N_n$ the number of brain regions, and $G_{ij}$ is the element of the gain matrix representing the distances of the sensors from the sources. The exponential function used in this study introduces a nonlinear relationship between source activity and sensor measurements, allowing for potentially capturing more complex dynamics. Assuming that the generated signal decays with the square of the distance from the source, the gain matrix is approximated by

$$G_{ij} = \sum_{k \in V_j} \frac{c A_k}{|\vec{x}_i^s - \vec{x}_k^v|^2}, \quad (3)$$

where $V_j$ is the set of all vertices on the triangulated surface of region $j$, $c$ is the scaling coefficient, $A_k$ is the surface area associated with vertex $k$, $\vec{x}_i^s$ is the position of the sensor $i$, and $\vec{x}_k^v$ is the position of the vertex $k$. Note that we have not taken into account the dependence of the
source-to-sensor decay on the orientation of the neuronal tissue, due to the lack of geometric information about the orientation in subcortical structures. Nevertheless, this consideration holds greater significance using neural fields compared to neural mass modeling [5, 10].

**State-space modeling.** State-space modeling [47–49] forms the fundamental basis of dynamical systems theory [50, 51] and control engineering [52, 53], to describe and analyze system dynamics over time, capturing interactions within data and actively manipulating and regulating system behavior to steer the system towards desired states or trajectories [54, 55]. Nonlinear state-space modeling further enhances the modeling of complex systems, for instance to capture seizure onset, progression, and offset, by incorporating nonlinear relationships and dynamics [13, 16, 17, 56].

In this study, the state-space representation of the VEP model is given by a system of non-linear differential equations as follows:

$$ \begin{cases} \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t), \boldsymbol{\theta}) + \mathbf{w}(t), & \mathbf{x}(0) = \mathbf{x}_{t_0} \\ \mathbf{y}(t) = \mathbf{h}(\mathbf{x}(t)) + \mathbf{v}(t) \end{cases} \qquad (4) $$

where $\mathbf{x}(t) \in \mathbb{R}^{N_n}$ is the $N_n$-dimensional vector of system states evolving over time, $\mathbf{x}_{t_0}$ is the initial state vector at time $t=0$, $\boldsymbol{\theta} \in \mathbb{R}^{N_p}$ contains all the unknown evolution parameters, $\mathbf{u}(t)$ stands for the external input, and $\mathbf{y}(t) \in \mathbb{R}^{N_s}$ denotes the measured data subject to the measurement error $\mathbf{v}(t)$. The process (dynamical) noise and the measurement noise denoted by $\mathbf{w}(t) \sim \mathcal{N}(0, \sigma^2)$ and $\mathbf{v}(t) \sim \mathcal{N}(0, \sigma'^2)$, respectively, are assumed to follow a Gaussian distribution with mean zero and variance $\sigma^2$ and $\sigma'^2$, respectively. Moreover, $\mathbf{f}(\cdot)$ is a vector function that describes the dynamical properties of the system i.e., summarizing the biophysical mechanisms underlying the temporal evolution of system states (here, govern by VEP model, Eq 1) and $\mathbf{h}(\cdot)$ represents a measurement function i.e., the instantaneous mapping from system states to observations (here, the gain matrix, Eqs 2, and 3).

**Parameter estimation problem.** In this study, it is assumed that the structure of the state-space model is known but that the associated model parameters $\boldsymbol{\theta} \in \{\eta_i, K\}$ with $i \in \{1, 2, \dots, N_n\}$ are unknown and need to be estimated from the available data. Using Desikan-Killiany parcellation with $N_n = 84$, $\boldsymbol{\theta} \in \mathbb{R}^{85}$. Although there are many strategies available for determining these parameters [17, 20, 21, 57–64], we focus on an approach based on global optimization. This method involves formulating an optimization problem that measures the difference between simulated data (produced by the model) and real observed data, aiming to adjust the model parameters for a better fit.

Given this context, to minimize the discrepancy between model predictions and observations, we employ the Root Mean Square Error (RMSE). A lower RMSE value indicates a better fit between model predictions and observed data, pointing to more accurate parameter estimation and enhanced model reliability.

Note that, due to the sparse placement of SEEG electrodes, the gain matrix is not of full rank (see S1 File), posing significant challenges for parameter estimation, particularly in accurately inferring the unknown combination of activity from neighboring brain regions near the sensor. These challenges include computational time, the reliability of estimating the epileptogenicity parameter, and the overall model inversion process. Consequently, the parameter estimation of dynamic models generally exhibits NP-hard complexity. Given the necessity to obtain a viable solution for the patient in a reasonable runtime, it is impractical to calculate for weeks to reach the optimal solution. Thus, in our global optimization approach, we have employed metaheuristics, local solvers, and HPC techniques, aiming to speedup the search as much as possible.
<table><caption>**Table 1. Optimization problems for various configurations of the VEP model.** The seizure propagation depends on the interplay between node dynamics (excitability) and network coupling (structure). The signals generated by VEP models on the brain region level are called source signals (see Eq 1). The measured signals from the electrodes are called sensor signals. To map the simulated sources from the brain regions to the sensors, the electromagnetic forward problem needs to be solved (see Eqs 2 and 3). Due to the sparse placement of electrodes, the gain (lead-field matrix) used for mapping the source to the sensors is not of full rank.</caption><thead><tr><td>problem</td><td>cost function</td><td>description</td></tr></thead><tbody><tr><td>1</td><td>deterministic</td><td>Forward simulation of VEP model with weak coupling at source level, i.e., no propagation, and the gain matrix G = I.</td></tr><tr><td>2</td><td>deterministic</td><td>Forward simulation of VEP model with weak coupling at sensor level, i.e., no propagation, and a low rank gain matrix.</td></tr><tr><td>3</td><td>deterministic</td><td>Forward simulation of VEP model with strong coupling at source level, i.e., with propagation, and the gain matrix G = I.</td></tr><tr><td>4</td><td>deterministic</td><td>Forward simulation of VEP model with strong coupling at sensor level, i.e., with propagation, and a low rank gain matrix.</td></tr><tr><td>5</td><td>stochastic</td><td>Forward simulation of SDE VEP model with large τ (stiff equations) at source level, with propagation, and the gain matrix G = I.</td></tr><tr><td>6</td><td>stochastic</td><td>Forward simulation of SDE VEP model with large τ (stiff equations) at sensor level, with propagation, and a low rank gain matrix.</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1011642.t001

Subsequently, we can rearrange various configurations of the VEP model for benchmarking based on the connectivity matrix of a patient. Table 1 presents 6 different optimization problems and a summary of their main characteristics. As we are assuming no noise in problems 1 to 4, they are treated as ODEs, and their parameter estimation problems are deterministic. However, in the case of problems 5 and 6, we handle the noise, so the SDEs result in our cost function in the calibration being stochastic.

**Cooperative optimization methods.** The domain of High-Performance Computing (HPC) offers essential methodologies and advanced technologies designed to significantly boost the efficiency of classical optimization algorithms, particularly within the context of contemporary many-core infrastructures. These methods and techniques play a pivotal role in specific applications where computation speed is intrinsically linked to the patient's well-being, underscoring the critical nature of rapid model calibration. In this context, messaging libraries like MPI emerge as key tools, providing an effective avenue for code parallelization. Concurrently, programming strategies such as the master-worker paradigm provide a streamlined and intuitive strategy for coordinating distributed optimization agents, further contributing to improve the overall efficacy and performance of numerical optimization methods.

Recently, leveraging the methodologies and technologies mentioned above, the Self-adaptive Cooperative enhanced Scatter Search (SaCeSS) [65] optimizer was presented as a competitive method to solve parameter estimation problems in large-scale nonlinear dynamic models. SaCeSS is a cooperative parallel (using an island-based model) implementation of an evolutionary algorithm, enhanced scatter search (eSS) [66, 67]. Independent eSS instances (workers) run in parallel in a cooperative fashion, and a master process manages communications between them. Each parallel eSS instance is a hybrid method, combining global search with calls to an efficient local search (LS) method (in this study, Dynamic Hill Climbing, DHC [68]).

During the execution of SaCeSS, if a worker explores a solution that improves the best-known one (considering all workers), this solution is asynchronously sent to the master, assessing its quality and determining whether it should be propagated to other MPI processes. This hybrid optimization strategy is enriched with several innovative mechanisms. These include: (i) the establishment of asynchronous collaboration amongst parallel processes,
<figure><img src="image_4.png" /><figcaption>**Fig 1.** General overview of SaCeSS. This method employs a parallel cooperative scheme based on a master-worker strategy. An example of solution propagation flow is illustrated: when a worker obtains a good solution, it is shared with the rest of the workers through the master. Additionally, the master implements an adaptation mechanism to replace the settings of those workers with poor performance (as registered in a scoreboard), helping to improve their performance.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g001

ensuring a seamless and efficient exchange of information; (ii) the integration of both coarse and fine-grained parallelism, providing a balanced and versatile computational structure; and (iii) the implementation of self-tuning strategies, which autonomously improve performance and adapt to varying conditions, ultimately enhancing the reliability and efficiency of the search process.

To provide a visual representation of the communication and adaptation mechanisms in SaCeSS, Fig 1 offers a detailed pictorial overview. In previous research, SaCeSS demonstrated good performance and robustness in addressing complex, large-scale model calibration problems within the domain of computational systems biology [69]. Leveraging these strengths, here we have tailored SaCeSS to tackle the specific optimization challenges associated with parameter estimation in VEP models.

**Scope and contributions.** In this study, we employed SaCeSS algorithm to calibrate various configurations of the VEP model, thereby demonstrating the potential of this solver for rapid parameter estimation in such contexts. As the VEP models serve as personalized brain
network models tailored to individual patients, they could lead to more informed clinical decision-making and optimize surgical strategies. For this reason, the primary objective is to ascertain whether our proposed approach can address the associated optimization problem within a reasonable computational time to avoid excessive waiting times for the patient.

We also evaluate the parallel efficiency and scalability of SaCeSS in the context of VEP models using two types of HPC infrastructures, a supercomputer and a PC workstation. Based on these results, we analyze the computational resources required to obtain fast and sufficiently accurate estimations, which is important to prevent extended delays in delivering medical results in a real-world application. Further, we also compare SaCeSS with other competitive parallel optimizers.

Additionally, since the calibration of large-scale models always contains some degree of uncertainty due to non-identifiabilities, we extend SaCeSS with a method of uncertainty quantification based on a recent ensemble approach [70, 71].

Finally, while SaCeSS is equipped with a self-adjusting mechanism for its settings, there remain certain global hyperparameters associated with its overall cooperation and adaptation functionalities. To enhance SaCeSS's robustness and efficiency even further, here we have used Optuna, a hyperparameter tuning tool [72]. This tool leverages Bayesian optimization techniques to determine the most effective global hyperparameters considering the specific nature of the VEP parameter estimation challenges.

## Results

In this section, our aim is to assess the performance of our proposed parameter estimation method in VEP models. To achieve this, we tested the SaCeSS method with the different benchmarks presented in Table 1. For each optimization problem, SaCeSS was run using the same computational resources: 12 parallel processors within a specific time threshold. Due to the non-deterministic nature of SaCeSS, it was necessary to perform multiple runs (10 times per problem) to ensure a comprehensive assessment. Recognizing the differing complexities of each benchmark and aiming to provide adequate time for convergence, we established the following stopping times (wall time limits) in the FT3 supercomputer: 3 hours for problem 1, 18 hours for problems 2–4, and 39 hours for problems 5 and 6.

Figs 2, 3 and 4 present the VEP configurations and estimations for different spatial maps of epileptogenicity, at source and sensors levels (see Table 1). For each problem, we have shown the whole-brain simulations (see Eq 1) at source level and corresponding SEEG simulations at sensor level (see Eq 2). Moreover, the trajectories in the phase-plane for different regions are illustrated. These results aim to showcase the capabilities of SaCeSS in capturing the true mechanism underlying seizure initiation and propagation from a dynamical systems theory perspective. We have utilized the confusion matrix, which was computed based on the inferred excitability $η_i$, to report the estimation accuracy for three node types: Epileptogenic Zone (EZ), Propagation Zone (PZ), and Healthy Zone (HZ).

Fig 2 shows the estimation result for the deterministic VEP model with weak coupling, i.e., no seizure recruitment from EZ to PZ. Fig 2A displays the simulated fast activities across brain regions at source level (cf., $x_i$ variables in Eq 1). In this scenario, two regions with high excitability $η_i$ are identified as EZ (shown in red), while three regions are designated as candidate regions with excitability close to the critical value for seizure propagation (PZ, shown in yellow). The remaining regions with low excitability are classified as healthy (HZ, shown in green). Fig 2B displays the true and estimated trajectories in the phase-plane. In the phase-planes depicted, the x- and z-nullclines are highlighted in blue, and the intersection of these nullclines indicates the fixed point of the system. A full circle represents a stable fixed point,
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Best solutions obtained by SaCeSS using 12 processors for the deterministic VEP model with weak coupling at source and sensor levels (problems 1 and 2).
(A) The simulated fast activities across brain regions at source level used for optimization. No propagation is observed due to weak coupling. (B) The true and estimated trajectories at phase-plane for three node types as Epileptogenic Zone (EZ, in red), Propagation Zone (PZ, in yellow), and Healthy Zone (HZ, in green). (C) The envelope of the simulated SEEG signals at the sensor level. (D) The confusion matrix indicates 100% accuracy in the estimation of three node types (EZ/PZ/HZ) based on SEEG signals.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g002

<figure><img src="image_6.png" /><figcaption>**Fig 3.** Best solutions obtained by SaCeSS using 12 processors for deterministic VEP model with strong coupling at source and sensor levels (problems 3 and 4). (A) The simulated fast activities across brain regions at source level used for optimization. Due to strong coupling, the seizure propagates from epileptogenic zone (in red) to other brain regions (in yellow). (B) The phase-plane displayed both the actual and predicted trajectories for three node categories as EZ, PZ, and HZ. (C) The simulated SEEG signals at the sensor level. (D) The confusion matrix that illustrates 100% accuracy for the estimation based on SEEG signals.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g003
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Best solutions obtained by SaCeSS using 12 processors for stochastic VEP model with stiff equations at source and sensor levels (problems 5 and 6). (A) The simulated fast activities across brain regions with large size lenght (due to large time scale separation $τ$) at source level used for optimization. Due to strong coupling, the seizure propagates from EZ (in red) to PZ (in yellow). (B) The true and estimated trajectories in phase-plane for EZ, PZ, and HZ. (C) The simulated SEEG signals at the sensor level. (D) The confusion matrix that illustrates that one of the PZ is mis-classified as HZ.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g004

while an empty circle represents an unstable fixed point. From this figure, we can observe that the trajectory of an HZ is attracted to the stable fixed point of the system (located on the left branch of the cubic x-nullcline), indicating that it does not trigger an epileptic seizure. For a PZ, the coupling and the value of excitability, which is close to the critical value of epileptogenicity, cause the z-nullcline to move downward. However, the coupling is not sufficient for propagation as the system's fixed point remains stable. For the EZ, the system exhibits an unstable fixed point due to the high value of excitability. In this regime, Epileptor possesses a limit cycle and the seizure triggers autonomously. We can observe that in all cases, the estimated trajectories closely follow the true trajectories due to the accurate estimation of excitability (see Table 2). Fig 2C shows the simulated SEEG signals at the sensor level given a low-rank sparse gain matrix. Fig 2D indicates a 100% accuracy in the estimation of three node types (EZ/PZ/HZ) based on SEEG signals. Overall, these results indicate that SaCeSS is able to accurately capture the system dynamics for different node types in the absence of noise in the VEP model.

Fig 3 shows the estimation result for the deterministic VEP model but with strong coupling, i.e., with seizure recruitment from EZ to PZ. Fig 3A displays the simulated fast activities across brain regions at the source level, where the seizure propagates to the candidate brain regions as PZ. In Fig 3B, the seizure recruitment to these regions is due to the network effects, as their equilibrium is close to the bifurcation value. Here, due to the sufficient coupling strength and the value of excitability which is close to the critical value of epileptogenicity, the z-nullcline moves down, causing a bifurcation thereby allowing the seizure to propagate. This indicates that seizure propagation depends on the interplay between the brain region's excitability (node dynamics), and the network coupling (parameter K). Fig 3C and 3D show the observed SEEG signals and the accuracy of the estimation. These results indicate that SaCeSS is able to
<table><caption>**Table 2.** Results for different parallel GO methods considering the set of VEP problems. Each solver was executed 10 times (different parallel jobs), with different stopping times depending on the problem. Resulted obtained using 12 parallel processors in the FT3 supercomputer. Here, fx stands for the cost function defined as the root mean square error between observed and generated data.</caption><thead><tr><th>problem</th><th>stopping time</th><th>best fx</th><th>mean fx ± std</th><th>problem</th><th>stopping time</th><th>best fx</th><th>mean fx ± std</th></tr></thead><tbody><tr><td colspan="4">SaCeSS</td><td colspan="4">PS-CMA-ES</td></tr><tr><td>1</td><td>3 h</td><td>5.30E-07</td><td>5.72E-07 ± 2.56E-08</td><td>1</td><td>3 h</td><td>1.70E-13</td><td>1.76E-13 ± 4.99E-15</td></tr><tr><td>2</td><td>18 h</td><td>2.20E-03</td><td>2.71E-03 ± 4.30E-04</td><td>2</td><td>18 h</td><td>1.70E+00</td><td>6.38E+00 ± 6.49E+00</td></tr><tr><td>3</td><td>18 h</td><td>4.93E-07</td><td>5.21E-07 ± 1.75E-08</td><td>3</td><td>18 h</td><td>1.59E-13</td><td>1.70E-13 ± 5.55E-15</td></tr><tr><td>4</td><td>18 h</td><td>3.39E-03</td><td>1.01E-02 ± 4.82E-03</td><td>4</td><td>18 h</td><td>2.82E+00</td><td>7.28E+00 ± 3.88E+00</td></tr><tr><td>5</td><td>39 h</td><td>7.35</td><td>7.74 ± 3.09E-01</td><td>5</td><td>39 h</td><td>22.48</td><td>24.46 ± 1.19</td></tr><tr><td>6</td><td>39 h</td><td>632.38</td><td>719.24 ± 130.47</td><td>6</td><td>39 h</td><td>1143.70</td><td>1302.89 ± 171.17</td></tr><tr><td colspan="4">asynPDE</td><td colspan="4">SaCeSS optuna</td></tr><tr><td>1</td><td>3 h</td><td>2.15E-05</td><td>8.68E-05 ± 1.01E-04</td><td>1</td><td>3 h</td><td>5.00E-07</td><td>5.47E-07 ± 3.31E-08</td></tr><tr><td>2</td><td>18 h</td><td>1.11E+00</td><td>1.35E+01 ± 1.37E+01</td><td>2</td><td>18 h</td><td>3.31E-04</td><td>7.32E-04 ± 3.15E-04</td></tr><tr><td>3</td><td>18 h</td><td>4.31E-13</td><td>5.99E-05 ± 9.18E-05</td><td>3</td><td>18 h</td><td>4.99E-07</td><td>5.15E-07 ± 1.17E-08</td></tr><tr><td>4</td><td>18 h</td><td>1.07</td><td>19.60 ± 9.47</td><td>4</td><td>18 h</td><td>4.90E-04</td><td>8.84E-04 ± 3.75E-04</td></tr><tr><td>5</td><td>39 h</td><td>27.01</td><td>30.95 ± 3.32</td><td>5</td><td>39 h</td><td>7.22</td><td>7.69 ± 0.25</td></tr><tr><td>6</td><td>39 h</td><td>961.30</td><td>1195.68 ± 109.48</td><td>6</td><td>39 h</td><td>535.88</td><td>792.44 ± 197.92</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1011642.t002

accurately capture the system dynamics for different node types with strong coupling in the VEP model.

Fig 4 shows the estimation result for the stochastic and stiff VEP model, i.e., with a long seizure envelope (see Fig 4A) due to a large time scale separation (parameter $\tau$ in Eq 1). In particular, we have considered a high value of noise dynamics, which is a zero-mean Gaussian noise with a standard deviation of 0.1. Therefore, given a known values of excitability and coupling, the seizure propagation is random and depends on the brain state dependency caused by the noise dynamics. Fig 4B shows that the system dynamics of regions corresponding to EZ and HZ are accurately estimated. However, the estimated trajectory corresponding to one of PZs dampens to a stable fixed point, despite the presence of a limit cycle in the observation. Fig 4C shows the simulated SEEG signals used for optimization, where one of the PZ regions is misclassified as HZ, as indicated by the confusion matrix (Fig 4D). This result indicates that while SaCeSS can effectively handle the very fast-changing components in the VEP model, accurately estimating the EZ/PZ from noisy signals can be challenging when using an error function such as RMSE for optimization.

Table 2 lists the best costs (RMSE) for each problem, alongside average and deviation values, where in problems from 1 to 4 (ODEs) approached zero. In contrast, the convergence of problems 5 and 6 (SDEs) is not close to zero as for the models based on ODEs, but the estimations are close to ground truth values.

We further assessed the scalability of SaCeSS on two platforms: the FinisTerrae III (FT3) supercomputer (details at https://www.cesga.es/en/), and a DELL Precision 5820 workstation equipped with 18 cores (Intel i9-10980XE at 3.00GHz). We conducted ten runs of each optimization problem using both the sequential eSS and the parallel SaCeSS solvers, employing 6, 12, and 24 processors—although only 6 and 12 cores were used in the case of the DELL Precision 5820. This exercise served to highlight the performance advantages of SaCeSS over eSS and offered valuable observations regarding the behavior of SaCeSS with varying numbers of parallel processors (workers).

The detailed results of this analysis are shown in Table 3, where, in general, all SaCeSS configurations achieved fitness functions close to zero, except for problems 5 and 6, where the
<table><caption>Table 3. Scalability of SaCeSS in an supercomputing infrastructure (FT3) and in a desktop workstation (DELL Precision 5820). SaCeSS was executed 10 times using different number of processors, varying the stopping time for each problem. In the case of one processor, the sequential enhanced scatter search (eSS) solver was used.</caption><thead><tr><th colspan="5">FT3 supercomputer</th><th colspan="5">DELL Precision 5820</th></tr><tr><th>problems</th><th>processors</th><th>stopping time</th><th>best fx</th><th>mean fx ± sd</th><th>problems</th><th>processors</th><th>stopping time</th><th>best fx</th><th>mean fx ± sd</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>1</td><td rowspan="4">3 h</td><td>5.96E-07</td><td>6.34E-07 ± 2.76E-08</td><td rowspan="4">1</td><td>1</td><td rowspan="4">3 h</td><td>1.11E-06</td><td>2.16E-06 ± 5.75E-07</td></tr><tr><td>6</td><td>5.37E-07</td><td>5.91E-07 ± 2.98E-08</td><td>6</td><td>6.42E-07</td><td>7.60E-07 ± 7.70E-08</td></tr><tr><td>12</td><td>5.30E-07</td><td>5.72E-07 ± 2.56E-08</td><td>12</td><td>5.88E-07</td><td>6.97E-07 ± 7.61E-08</td></tr><tr><td>24</td><td>5.10E-07</td><td>5.47E-07 ± 1.74E-08</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">2</td><td>1</td><td rowspan="4">18 h</td><td>8.20E-03</td><td>7.63E-02 ± 2.02E-01</td><td rowspan="4">2</td><td>1</td><td rowspan="4">3 h</td><td>3.31E-02</td><td>5.92E+00 ± 7.63E+00</td></tr><tr><td>6</td><td>2.78E-03</td><td>3.28E-03 ± 3.78E-04</td><td>6</td><td>6.35E-03</td><td>1.52E+00 ± 2.30E+00</td></tr><tr><td>12</td><td>2.20E-03</td><td>2.71E-03 ± 4.30E-04</td><td>12</td><td>3.54E-03</td><td>5.79E-03 ± 2.51E-03</td></tr><tr><td>24</td><td>1.31E-03</td><td>1.77E-03 ± 4.22E-04</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">3</td><td>1</td><td rowspan="4">18 h</td><td>5.12E-07</td><td>5.63E-07 ± 3.01E-08</td><td rowspan="4">3</td><td>1</td><td rowspan="4">3 h</td><td>5.42E-07</td><td>7.55E-07 ± 1.37E-07</td></tr><tr><td>6</td><td>5.02E-07</td><td>5.34E-07 ± 1.76E-08</td><td>6</td><td>5.28E-07</td><td>5.71E-07 ± 2.48E-08</td></tr><tr><td>12</td><td>4.93E-07</td><td>5.21E-07 ± 1.75E-08</td><td>12</td><td>5.08E-07</td><td>5.66E-07 ± 3.40E-08</td></tr><tr><td>24</td><td>4.91E-07</td><td>5.14E-07 ± 1.15E-08</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">4</td><td>1</td><td rowspan="4">18 h</td><td>6.86E-03</td><td>1.34E-02 ± 2.89E-03</td><td rowspan="4">4</td><td>1</td><td rowspan="4">3 h</td><td>2.07E-02</td><td>2.25E+01 ± 6.08E+01</td></tr><tr><td>6</td><td>4.83E-03</td><td>1.15E-02 ± 6.45E-03</td><td>6</td><td>1.74E-02</td><td>3.49E+00 ± 7.55E+00</td></tr><tr><td>12</td><td>3.39E-03</td><td>1.01E-02 ± 4.82E-03</td><td>12</td><td>4.51E-03</td><td>1.40E+00 ± 3.13E+00</td></tr><tr><td>24</td><td>1.49E-03</td><td>2.03E-03 ± 3.83E-04</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">5</td><td>1</td><td rowspan="4">39 h</td><td>8.00</td><td>8.38 ± 0.23</td><td rowspan="4">5</td><td>1</td><td rowspan="4">3 h</td><td>11.99</td><td>13.26 ± 0.96</td></tr><tr><td>6</td><td>7.51</td><td>7.89 ± 0.28</td><td>6</td><td>10.10</td><td>10.71 ± 0.43</td></tr><tr><td>12</td><td>7.34</td><td>7.73 ± 0.30</td><td>12</td><td>9.71</td><td>10.07 ± 0.21</td></tr><tr><td>24</td><td>6.93</td><td>7.45 ± 0.32</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="4">6</td><td>1</td><td rowspan="4">39 h</td><td>691.45</td><td>975.04 ± 177.04</td><td rowspan="4">6</td><td>1</td><td rowspan="4">3 h</td><td>918.34</td><td>1086.80 ± 82.48</td></tr><tr><td>6</td><td>641.97</td><td>777.83 ± 171.07</td><td>6</td><td>782.74</td><td>990.37 ± 127.55</td></tr><tr><td>12</td><td>632.38</td><td>719.24 ± 130.47</td><td>12</td><td>730.25</td><td>948.03 ± 142.55</td></tr><tr><td>24</td><td>593.41</td><td>638.19 ± 29.35</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1011642.t003

RMSE metric is unable to capture the noisy nature of the observation at the sensor level. The results indicate that increasing the number of processors (workers) leads to an enhancement in the accuracy of solutions. This improvement is also mirrored in the reduction of mean and deviation values. Fig 5 provides additional details: the boxplots show the spread of the best solutions obtained, and the convergence curves (evolution of objective function versus time considering the best run of each method) illustrate the speedup gains with more processors. The bootstrapping on different runs complements this observation, indicating that the distribution of solutions has reduced its variability and is pushed towards better results.

On the other hand, we also evaluated the behavior of SaCeSS when scaling the dimension of the problem. To investigate this, we generated instances of varying sizes using the Desikan-Killiany parcellation with 42, 84 regions, the VEP atlas with 162 regions, and Schaefer atlas with 400 regions, for problem 1, which is the most tractable for this experiment with all parameters identifiable. We then proceeded to solve these instances using SaCeSS with 12 processors and a time threshold of 3 hours. In Fig 6, the convergence curves and solution boxplots indicate that although the dimension increases and convergence takes longer to reach RMSE values close to zero, SaCeSS performs reasonably well, and no run gets stuck at an unusual value.

Additionally, we benchmarked SaCeSS against two competitive parallel global optimizers which have previously shown good performance in various computational biology problems: the Particle Swarm CMA Evolution Strategy (PS-CMA-ES) [73] and the asynchronous parallel
<figure><img src="image_8.png" /><figcaption>Problem 1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Problem 2</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Problem 3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Problem 4</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Problem 5</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Problem 6</figcaption></figure>

**Fig 5. Scalability analysis of SaCeSS on the Finisterrae III supercomputer.** Each convergence curve represents the best run for the sequential method (eSS) and for SaCeSS using 6, 12, and 24 processors. The boxplots illustrate the spread in the solutions obtained in repeated runs (using different number of parallel workers for the case of SaCeSS).

https://doi.org/10.1371/journal.pcbi.1011642.g005

Differential Evolution (asynPDE) [74]. Maintaining the same computational testbed with 12 parallel processors per run and the previously mentioned time thresholds, we analyzed for these solvers the quality of the best solutions and the associated means and deviations. Table 2 presents a comprehensive summary of results, and Fig 7 shows the convergence curves and final solution spread (with boxplots). These results unequivocally demonstrate the enhanced performance of our proposed method over both asynPDE and PS-CMA-ES in various aspects of solving VEP problems. Our method exhibits a more rapid convergence rate, better robustness in the estimations (as indicated by a lower dispersion in results across multiple runs), and the ability to find superior solutions. The sole exceptions are observed when addressing the deterministic VEP at the source level (problems 1 and 3), where the estimations are much easier due to the full identifiability of all parameters involved. In these particular problems, PS-CMA-ES yields results very close to zero (in the order of $1E-13$), outperforming SaCeSS, which provides results in the order of $1E-7$. While SaCeSS converges more swiftly, it tends to level off at these values. Nonetheless, it is important to note that these fits are exceptionally accurate and, from a practical standpoint, they are virtually indistinguishable.
<figure><img src="image_14.png" /><figcaption>**Fig 6. Scalability analysis of SaCeSS with different numbers of brain regions.** The four groups of convergence curves (differentiated by color) represent runs for SaCeSS using 12 processors solving problem 1, VEP with weak coupling at the source level, having different number of brain regions (NN): 42, 84, 162, and 400. The boxplots illustrate the spread in the solutions obtained in repeated runs.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g006

We also explored the use of a Bayesian optimization framework (via Optuna [72]) to further improve the performance of SaCeSS in the VEP benchmarks. The objective was to fine-tune specific configuration settings related to the cooperation mechanism. This approach required the formulation of a new mixed-integer optimization problem with six decision variables associated with the configuration options to be tuned. The objective was to minimize a cost function defined by the geometric mean of cost values from five independent SaCeSS runs. Each of these five runs used 12 processors with stopping criteria only based on time (30 seconds). We used this threshold because it is the time required to reach quality solutions in problem 1. In other words, we used the Bayesian optimization scheme provided by Optuna to obtain the best SaCeSS hyperparameters for problem 1. This process required 13 hours to complete and was repeated 10 times due to the stochastic nature of Bayesian optimization. We only performed this Bayesian tuning for problem 1, the easiest case study. Detailed results are given in S1 File. The computational cost of repeating this fine-tuning for the other problems would be very significant.

For this reason, we applied the best settings obtained during the calibration of problem 1 to the others, assuming that the Bayesian fine-tuning would generalize well. The results are reported in Table 2, where the fine-tuned SaCeSS outperformed the original SaCeSS configuration in terms of both the best solution and its mean. Fig 8 provides a side-by-side comparison, visually illustrating the differences in solution spread and convergence curves (across all runs) between the standard and the fine-tuned SaCeSS. The latter converged faster and to better solutions in problems 1–4, reduced the spread in the results. However, the results for problems 5 and 6 were not statistically different, indicating that the SaCeSS default self-adaptive mechanisms are capable of handling the more challenging problems.

Regarding uncertainty quantification, we used an approach based on ensemble modeling [70, 71]. This method takes samples from the SaCeSS optimization in the vicinity of the global solution, building an ensemble of calibrations which represent equally well the data.
<figure><img src="image_15.png" /><figcaption>**Fig 7.** Comparison of SaCeSS with other parallel optimizers. Each convergence curve represents the best run for SaCeSS, asynPDE, and PS-CMA-ES, all using 12 processors. The boxplots illustrate the dispersion in the solutions obtained with the aforementioned solvers. Each colored spot within the boxplots represents the solution cost obtained by a single run for each solver.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g007

During the optimizations, SaCeSS stores all the explored solutions of certain quality (defined by a threshold). After convergence, our workflow selects a representative subsample of the stored solutions, normalizes it, and then illustrates the variability in the estimated parameters for each problem using parallel coordinate graphs. These results are shown in Fig 9, wherein the y-axis, the time-varying distribution percentiles are shown as shaded red bands around a central black line (the median). The x-axis represents the parameters' index in the VEP models: the level of epileptogenicity $\eta_i$ at each brain region and the global coupling K (for the label names and indices of the sub-divided brain regions, see S1 File). Due to the sparsity of the gain matrix, we observe a higher level of uncertainty in the estimation at the sensor level (problems 2, 4, and 6) compared to the source level (problems 1, 3, and 5).

## Discussion

The accurate parameter estimations in different configurations of the VEP model, as shown in Table 2 and Figs 2, 3 and 4 highlight the capability of SaCeSS to solve high-dimensional problems within a reasonable computational time, using just 12 processors (additional details can be found at https://doi.org/10.5281/zenodo.10057788). The convergence curves further substantiate this observation: in different deterministic problems, SaCeSS exhibit fast convergence while obtaining solutions close to the global optimum. The use of SaCeSS with 12 processors and a time threshold of 3 hours demonstrates the algorithm's convergence in deterministic
<figure><img src="image_16.png" /><figcaption>**Fig 8.** Comparison between the original and the fine-tuned SaCeSS. Convergence curves for original (orange) and fine-tuned SaCeSS (blue), plus boxplots illustrating the solutions spread.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g008

problems, even for a large number of regions (see Fig 6), confirming the method's scalability. However, in stochastic problems, SaCeSS presents a less efficient behavior due to the complexity introduced by dynamical noise. In particular, the solution for problem 6 did not estimate the phase-plane trajectories in the PZ region adequately. Nevertheless, considering a slightly modified cost function for problem 6 (details given in S1 File), SaCeSS was able to obtain better fitness than the ground truth model. Technically this corresponds to a slight overfitting, thus explaining why we could not recover the nominal values of the parameters. However, we also suspect that there are identifiability issues caused by lack of enough information in the data. Thus, in order to enhance efficiency when addressing stochastic models, a promising strategy could be to study the integration of SaCeSS with Bayesian optimization methods, known for their proficiency in such challenges, at least for problems of moderate number of parameters.

Moreover, our method demonstrates versatility across different computational infrastructures. The original SaCeSS was designed for high-performance computing environments because of its parallel nature and scalability, i.e. adding parallel workers to the search increases the probability of finding a better solution faster. Nevertheless, SaCeSS also delivers steadfast performance on smaller-scale systems, such as desktop workstations equipped with a moderate
<figure><img src="image_17.png" /><figcaption>**Fig 9.** Uncertainty quantification of the parameters estimated using SaCeSS: Dispersion as parallel coordinates plots. On the y-axis, the time-varying distribution percentiles are manifested as shaded red bands encircling a central black line, which signifies the median. The x-axis represents the parameters in the VEP models: the level of epileptogenicity ηᵢ at each brain region and the global coupling K.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011642.g009
quantity of cores and memory, even when operating under stringent time limitations (see Table 3, “DELL Precision 5820” results). This flexibility ensures accessibility for a broad user base, including those with limited experience in high-performance computing, allowing them to fully exploit the advantages offered by our methodology.

In the comparison of SaCeSS with two other leading parallel solvers (asynPDE and PS-CMA-ES), as showcased in Table 2, our approach demonstrates a consistent advantage over these competitors across multiple criteria in addressing VEP problems. These results highlight the competitive edge of cooperative parallel methods like SaCeSS against other strategies utilizing more conventional parallelization techniques. Moreover, considering the scarcity of parallel global optimization methods, and the complexities entailed in handling estimation problems involving differential equations, this comparative analysis emerges as a significant and informative contribution in itself.

Finally, the hyperparameter fine-tuning of SaCeSS using Bayesian optimization (with the Optuna solver) has resulted in significant improvements, particularly in the VEP problems with a deterministic cost function. Despite the time-consuming nature of evaluating the associated cost function, Optuna can provide a good solution without requiring excessive iterations. The fine-tuning revealed that four of the hyperparameters are similar to their default values, which were established based on our previous experience with SaCeSS optimizations. However, the Optuna results found larger values for two specific settings: the maximum number of evaluations for local solvers and the maximum number of evaluations before reconfiguring inactive workers. These larger values result in a more diverse search strategy, allowing local and global solvers more opportunities to identify and propagate effective solutions. We regard this as a first step exploring a promising field: the fine-tuning of cooperative parallel methods for challenging global optimization problems. In this particular context, as opposed to sequential algorithms, the settings to be fine-tuned are related to “social” features of multi-agent systems, such as communications among different workers. Therefore, several open questions remain for future consideration, including the calibration of stochastic problems given their high computational cost, and enhancing this procedure through a sensitivity analysis to determine the most impactful settings for inclusion or exclusion during fine-tuning.

In summary, driven by the challenge of estimating parameters in a series of progressively complex whole-brain network models related to epilepsy, this study started with the objective of selecting a resilient and efficient global optimization solver. As a second objective, we wanted a solver capable of leveraging parallel computing across various platforms, ranging from personal desktop workstations to supercomputers.

Our investigation revealed the robustness and efficiency of the SaCeSS parallel solver in addressing the complex parameter estimation problems associated with VEP models. This algorithm outperformed other parallel solvers in most VEP benchmark problems and exhibited enhanced performance when scaling the number of parallel processors. Additionally, using a Bayesian framework for hyperparameter tuning, we were able to further improve its performance. We also extended SaCeSS with a mechanism to store and post-process the parameter space sampled in the vicinity of the global solution and, after convergence, use the sampling to perform approximate uncertainty quantification. This extension allowed a computationally efficient assessment of the variability in the estimated parameters. Overall, these findings collectively highlight the potential of our global optimization approach in contributing to more informed clinical decision-making through fast and accurate parameter estimation in VEP models. However, this needs to be investigated using real empirical data in future studies.
**Supporting information**

S1 File. Supporting material document.

(PDF)

**Acknowledgments**

The authors acknowledge CESGA (Centro de Supercomputación de Galicia) for providing access to its FinisTerrae III supercomputer.

**Author Contributions**

**Conceptualization:** Meysam Hashemi, Viktor K. Jirsa, Julio R. Banga.

**Data curation:** Meysam Hashemi.

**Formal analysis:** Meysam Hashemi.

**Funding acquisition:** Viktor K. Jirsa, Julio R. Banga.

**Methodology:** David R. Penas, Meysam Hashemi, Julio R. Banga.

**Software:** David R. Penas, Meysam Hashemi.

**Supervision:** Viktor K. Jirsa, Julio R. Banga.

**Visualization:** David R. Penas, Meysam Hashemi.

**Writing – original draft:** David R. Penas, Meysam Hashemi.

**Writing – review & editing:** David R. Penas, Meysam Hashemi, Viktor K. Jirsa, Julio R. Banga.

**References**

1. Talairach J, Bancaud J. Lesion, irritative zone and epileptogenic focus. Confinia neurologica. 1966; 27 (1-3):91-94. https://doi.org/10.1159/000103937 PMID: 5334025

2. Bancaud J, Angelergues R, Bernouilli C, Bonis A, Bordas-Ferrer M, Bresson M, et al. Functional stereotaxic exploration (SEEG) of epilepsy. Electroencephalography and clinical neurophysiology. 1970; 28 (1):85-86. PMID: 4188481

3. Rosenow F, Lüders H. Presurgical evaluation of epilepsy. Brain. 2001; 124(9):1683-1700. https://doi.org/10.1093/brain/124.9.1683 PMID: 11522572

4. Makhalova J, Medina Villalon S, Wang H, Giusiano B, Woodman M, Bénar C, et al. Virtual epileptic patient brain modeling: Relationships with seizure onset and surgical outcome. Epilepsia. 2022; 63 (8):1942-1955. https://doi.org/10.1111/epi.17310 PMID: 35604575

5. Jirsa V, Wang H, Triebkorn P, Hashemi M, Jha J, Gonzalez-Martinez J, et al. Personalised virtual brain models in epilepsy. The Lancet Neurology. 2023; 22(5):443-454. https://doi.org/10.1016/S1474-4422(23)00008-X PMID: 36972720

6. Bartolomei F, Lagarde S, Wendling F, McGonigal A, Jirsa V, Guye M, et al. Defining epileptogenic networks: contribution of SEEG and signal analysis. Epilepsia. 2017; 58(7):1131-1147. https://doi.org/10.1111/epi.13791 PMID: 28543030

7. Proix T, Bartolomei F, Guye M, Jirsa VK. Individual brain structure and modelling predict seizure propagation. Brain. 2017; 140(3):641-654. https://doi.org/10.1093/brain/awx004 PMID: 28364550

8. Englot DJ, Nagarajan SS, Imber BS, Raygor KP, Honma SM, Mizuiri D, et al. Epileptogenic zone localization using magnetoencephalography predicts seizure freedom in epilepsy surgery. Epilepsia. 2015; 56(6):949-958. https://doi.org/10.1111/epi.13002 PMID: 25921215

9. Lopes M, Goodfellow M, Terry J. A model-based assessment of the seizure onset zone predictive power to inform the epileptogenic zone. Front Comput Neurosci. 2019; 13(25). https://doi.org/10.3389/fncom.2019.00025 PMID: 31105545
10. Wang HE, Woodman M, Triebkorn P, Lemarechal JD, Jha J, Dollomaja B, et al. Delineating epilepto- genic networks using brain imaging data and personalized modeling in drug-resistant epilepsy. Science Translational Medicine. 2023; 15(680):eabp8982. https://doi.org/10.1126/scitranslmed.abp8982 PMID: 36696482

11. Jobst BC, Cascino GD. Resective epilepsy surgery for drug-resistant focal epilepsy: a review. Jama. 2015; 313(3):285-293. https://doi.org/10.1001/jama.2014.17426 PMID: 25602999

12. Kuhlmann L, Lehnertz K, Richardson MP, Schelter B, Zaveri HP. Seizure prediction-ready for a new era. Nature Reviews Neurology. 2018; 14(10):618-630. https://doi.org/10.1038/s41582-018-0055-2 PMID: 30131521

13. Hashemi M, Vattikonda AN, Sip V, Diaz-Pier S, Peyser A, Wang H, et al. On the influence of prior information evaluated by fully Bayesian criteria in a personalized whole-brain model of epilepsy spread. PLoS computational biology. 2021; 17(7):e1009129. https://doi.org/10.1371/journal.pcbi.1009129 PMID: 34260596

14. Cao M, Galvis D, Vogrin SJ, Woods WP, Vogrin S, Wang F, et al. Virtual intracranial EEG signals reconstructed from MEG with potential for epilepsy surgery. Nature Communications. 2022; 13(1):994. https://doi.org/10.1038/s41467-022-28640-x PMID: 35194035

15. Sohrabpour A, Cai Z, Ye S, Brinkmann B, Worrell G, He B. Noninvasive electromagnetic source imaging of spatiotemporally distributed epileptogenic brain sources. Nature communications. 2020; 11 (1):1946. https://doi.org/10.1038/s41467-020-15781-0 PMID: 32327635

16. Jirsa VK, Proix T, Perdikis D, Woodman MM, Wang H, Gonzalez-Martinez J, et al. The Virtual Epileptic Patient: individualized whole-brain models of epilepsy spread. Neuroimage. 2017; 145:377-388. https://doi.org/10.1016/j.neuroimage.2016.04.049 PMID: 27477535

17. Hashemi M, Vattikonda A, Sip V, Guye M, Bartolomei F, Woodman MM, et al. The Bayesian Virtual Epileptic Patient: A probabilistic framework designed to infer the spatial map of epileptogenicity in a personalized large-scale brain model of epilepsy spread. NeuroImage. 2020; 217:116839. https://doi.org/10.1016/j.neuroimage.2020.116839 PMID: 32387625

18. Kini LG, Bernabei JM, Mikhail F, Hadar P, Shah P, Khambhati AN, et al. Virtual resection predicts surgical outcome for drug-resistant epilepsy. Brain. 2019; 142(12):3892-3905. https://doi.org/10.1093/brain/awz303 PMID: 31599323

19. Sip V, Hashemi M, Vattikonda AN, Woodman MM, Wang H, Scholly J, et al. Data-driven method to infer the seizure propagation patterns in an epileptic brain from intracranial electroencephalography. PLoS computational biology. 2021; 17(2):e1008689. https://doi.org/10.1371/journal.pcbi.1008689 PMID: 33596194

20. Vattikonda AN, Hashemi M, Sip V, Woodman MM, Bartolomei F, Jirsa VK. Identifying spatio-temporal seizure propagation patterns in epilepsy using Bayesian inference. Communications biology. 2021; 4 (1):1244. https://doi.org/10.1038/s42003-021-02751-5 PMID: 34725441

21. Jha J, Hashemi M, Vattikonda AN, Wang H, Jirsa V. Fully Bayesian estimation of virtual brain parameters with self-tuning Hamiltonian Monte Carlo. Machine Learning: Science and Technology. 2022; 3 (3):035016.

22. Hashemi M, Vattikonda AN, Jha J, Sip V, Woodman MM, Bartolomei F, et al. Amortized Bayesian inference on generative dynamical network models of epilepsy using deep neural density estimators. Neural Networks. 2023; 163:178-194. https://doi.org/10.1016/j.neunet.2023.03.040 PMID: 37060871

23. Schirner M, Rothmeier S, Jirsa VK, McIntosh AR, Ritter P. An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data. NeuroImage. 2015; 117:343-357. https://doi.org/10.1016/j.neuroimage.2015.03.055 PMID: 25837600

24. Fischl B. FreeSurfer. NeuroImage. 2012; 62(2):774-781. https://doi.org/10.1016/j.neuroimage.2012.01.021 PMID: 22248573

25. Jenkinson M, Bannister P, Brady M, Smith S. Improved Optimization for the Robust and Accurate Linear Registration and Motion Correction of Brain Images. NeuroImage. 2002; 17(2):825-841. https://doi.org/10.1006/nimg.2002.1132 PMID: 12377157

26. Tournier JD, Calamante F, Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007; 35 (4):1459-1472. https://doi.org/10.1016/j.neuroimage.2007.02.016 PMID: 17379540

27. Tournier JD, Calamante F, Connelly A. Determination of the appropriate value and number of gradient directions for high-angular-resolution diffusion-weighted imaging. NMR in Biomedicine. 2013; 26 (12):1775-1786. https://doi.org/10.1002/nbm.3017 PMID: 24038308

28. Tournier JD, Calamante F, Connelly A. Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions. In: Proceedings of the international society for magnetic resonance in medicine. vol. 18; 2010. p. 1670.
29. Desikan RS, Ségonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. *NeuroImage*. 2006; 31(3):968-980. https://doi.org/10.1016/j.neuroimage.2006.01.021 PMID: 16530430

30. Sanz-Leon P, Knock SA, Spiegler A, Jirsa VK. Mathematical framework for large-scale brain network modeling in The Virtual Brain. *Neuroimage*. 2015; 111:385-430. https://doi.org/10.1016/j.neuroimage.2015.01.002 PMID: 25592995

31. Jirsa VK, Haken H. Field Theory of Electromagnetic Brain Activity. *Phys Rev Lett*. 1996; 77(5):960-963. https://doi.org/10.1103/PhysRevLett.77.960 PMID: 10062950

32. David O, Friston KJ. A neural mass model for MEG/EEG: coupling and neuronal dynamics. *NeuroImage*. 2003; 20(3):1743-1755. https://doi.org/10.1016/j.neuroimage.2003.07.015 PMID: 14642484

33. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston K. The dynamic brain: from spiking neurons to neural masses and cortical fields. *PloS Comp Biol*. 2008; 4(8). https://doi.org/10.1371/journal.pcbi.1000092 PMID: 18769680

34. Cook BJ, Peterson AD, Woldman W, Terry JR. Neural Field Models: A mathematical overview and unifying framework. *Mathematical Neuroscience and Applications*. 2022; 2. https://doi.org/10.46298/mna.7284

35. Hashemi M, Hutt A, Sleigh J. Anesthetic action on extra-synaptic receptors: effects in neural population models of EEG activity. *Frontiers in Systems Neuroscience*. 2014; 8:232. https://doi.org/10.3389/fnsys.2014.00232 PMID: 25540612

36. Courtiol J, Guye M, Bartolomei F, Petkoski S, Jirsa VK. Dynamical Mechanisms of Interictal Resting-State Functional Connectivity in Epilepsy. *Journal of Neuroscience*. 2020; 40(29):5572-5588. https://doi.org/10.1523/JNEUROSCI.0905-19.2020 PMID: 32513827

37. Lavanga M, Stumme J, Yalcinkaya BH, Fousek J, Jockwitz C, Sheheitli H, et al. The virtual aging brain: a model-driven explanation for cognitive decline in older subjects. *bioRxiv*. 2022;.

38. Yalcinkaya BH, Ziaeemehr A, Fousek J, Hashemi M, Lavanga M, Solodkin A, et al. Personalized virtual brains of Alzheimer's Disease link dynamical biomarkers of fMRI with increased local excitability. *medRxiv*. 2023; p. 2023-01.

39. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. *Brain*. 2014; 137(8):2210-2230. https://doi.org/10.1093/brain/awu133 PMID: 24919973

40. El Houssaini K, Bernard C, Jirsa VK. The epileptor model: a systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block. *Eneuro*. 2020; 7(2). https://doi.org/10.1523/ENEURO.0485-18.2019 PMID: 32066612

41. Saggio ML, Crisp D, Scott JM, Karoly P, Kuhlmann L, Nakatani M, et al. A taxonomy of seizure dynamo-types. *Elife*. 2020; 9:e55632. https://doi.org/10.7554/eLife.55632 PMID: 32691734

42. Haken H. *Synergetics*. Physics Bulletin. 1977; 28(9):412. https://doi.org/10.1088/0031-9112/28/9/027

43. Jirsa VK, Friedrich R, Haken H, Kelso JS. A theoretical model of phase transitions in the human brain. *Biological cybernetics*. 1994; 71:27-35. https://doi.org/10.1007/BF00198909 PMID: 8054384

44. Hashemi M, Valizadeh A, Azizi Y. Effect of duration of synaptic activity on spike rate of a Hodgkin-Huxley neuron with delayed feedback. *Physical Review E*. 2012; 85(2):021917. https://doi.org/10.1103/PhysRevE.85.021917 PMID: 22463254

45. Jafarian A, Zeidman P, Wykes RC, Walker M, Friston KJ. Adiabatic dynamic causal modelling. *NeuroImage*. 2021; 238:118243. https://doi.org/10.1016/j.neuroimage.2021.118243 PMID: 34116151

46. Olmi S, Petkoski S, Guye M, Bartolomei F, Jirsa V. Controlling seizure propagation in large-scale brain networks. *PLoS computational biology*. 2019; 15(2):e1006805. https://doi.org/10.1371/journal.pcbi.1006805 PMID: 30802239

47. Durbin J, Koopman SJ. Time series analysis by state space methods. vol. 38. OUP Oxford; 2012.

48. Box GE, Jenkins GM, Reinsel GC, Ljung GM. Time series analysis: forecasting and control. John Wiley & Sons; 2015.

49. Hashemi M, Hutt A, Buhry L, Sleigh J. Optimal model parameter estimation from EEG power spectrum features observed during general anesthesia. *Neuroinformatics*. 2018; 16:231-251. https://doi.org/10.1007/s12021-018-9369-x PMID: 29516302

50. Ionides EL, Bretó C, King AA. Inference for nonlinear dynamical systems. *Proceedings of the National Academy of Sciences*. 2006; 103(49):18438-18443. https://doi.org/10.1073/pnas.0603181103 PMID: 17121996

51. Izhikevich EM. Dynamical systems in neuroscience. MIT press; 2007.

52. Trentelman HL, Stoorvogel AA, Hautus M. Control theory for linear systems. Springer Science & Business Media; 2012.
53. Medaglia JD, Pasqualetti F, Hamilton RH, Thompson-Schill SL, Bassett DS. Brain and cognitive reserve: Translation via network control theory. Neuroscience & Biobehavioral Reviews. 2017; 75:53-64. https://doi.org/10.1016/j.neubiorev.2017.01.016 PMID: 28104411

54. Carver CS, Scheier MF. Control theory: A useful conceptual framework for personality-social, clinical, and health psychology. Psychological bulletin. 1982; 92(1):111. https://doi.org/10.1037/0033-2909.92.1.111 PMID: 7134324

55. Scheid BH, Ashourvan A, Stiso J, Davis KA, Mikhail F, Pasqualetti F, et al. Time-evolving controllability of effective connectivity networks during seizure progression. Proceedings of the National Academy of Sciences. 2021; 118(5):e2006436118. https://doi.org/10.1073/pnas.2006436118 PMID: 33495341

56. Proix T, Jirsa VK, Bartolomei F, Guye M, Truccolo W. Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy. Nature communications. 2018; 9(1):1088. https://doi.org/10.1038/s41467-018-02973-y PMID: 29540685

57. Ghahramani Z, Hinton GE. Variational learning for switching state-space models. Neural computation. 2000; 12(4):831-864. https://doi.org/10.1162/089976600300015619 PMID: 10770834

58. Turner R, Deisenroth M, Rasmussen C. State-space inference and learning with Gaussian processes. In: Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. JMLR Workshop and Conference Proceedings; 2010. p. 868-875.

59. Frigola R, Chen Y, Rasmussen CE. Variational Gaussian process state-space models. Advances in neural information processing systems. 2014;27.

60. Archer E, Park IM, Buesing L, Cunningham J, Paninski L. Black box variational inference for state space models. arXiv preprint arXiv:151107367. 2015;

61. Brunton SL, Proctor JL, Kutz JN. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences. 2016; 113(15):3932-3937. https://doi.org/10.1073/pnas.1517384113 PMID: 27035946

62. Nassar J, Linderman SW, Bugallo M, Park IM. Tree-structured recurrent switching linear dynamical systems for multi-scale modeling. arXiv preprint arXiv:181112386. 2018;

63. Koppe G, Toutounji H, Kirsch P, Lis S, Durstewitz D. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI. PLoS computational biology. 2019; 15(8):e1007263. https://doi.org/10.1371/journal.pcbi.1007263 PMID: 31433810

64. Sip V, Hashemi M, Dickscheid T, Amunts K, Petkoski S, Jirsa V. Characterization of regional differences in resting-state fMRI with a data-driven network model of brain dynamics. Science Advances. 2023; 9(11):eabq7547. https://doi.org/10.1126/sciadv.abq7547 PMID: 36930710

65. Penas DR, González P, Egea JA, Doallo R, Banga JR. Parameter estimation in large-scale systems biology models: a parallel and self-adaptive cooperative strategy. BMC bioinformatics. 2017; 18:1-24. https://doi.org/10.1186/s12859-016-1452-4 PMID: 28109249

66. Rodriguez-Fernandez M, Egea JA, Banga JR. Novel Metaheuristic for Parameter Estimation in Nonlinear Dynamic Biological Systems. BMC Bioinformatics. 2006; 7(1):483. https://doi.org/10.1186/1471-2105-7-483 PMID: 17081289

67. Egea JA, Martí R, Banga JR. An evolutionary method for complex-process optimization. Computers & Operations Research. 2010; 37(2):315-324. https://doi.org/10.1016/j.cor.2009.05.003

68. de la Maza M, Yuret D. Dynamic hill climbing. AI Expert. 1994; 9(3):26-31.

69. Penas DR, Henriques D, González P, Doallo R, Saez-Rodriguez J, Banga JR. A parallel metaheuristic for large mixed-integer dynamic optimization problems, with applications in computational biology. PLOS ONE. 2017; 12(8):1-32. https://doi.org/10.1371/journal.pone.0182186 PMID: 28813442

70. Villaverde AF, Bongard S, Mauch K, Müller D, Balsa-Canto E, Schmid J, et al. A consensus approach for estimating the predictive accuracy of dynamic models in biology. Computer methods and programs in biomedicine. 2015; 119(1):17-28. https://doi.org/10.1016/j.cmpb.2015.02.001 PMID: 25716416

71. Villaverde AF, Raimúndez E, Hasenauer J, Banga JR. Assessment of prediction uncertainty quantification methods in systems biology. IEEE/ACM transactions on computational biology and bioinformatics. 2022;

72. Akiba T, Sano S, Yanase T, Ohta T, Koyama M. Optuna: A next-generation hyperparameter optimization framework. In: Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining; 2019. p. 2623-2631.

73. Müller CL, Baumgartner B, Ofenbeck G, Schrader B, Sbalzarini IF. pCMALib: a parallel fortran 90 library for the evolution strategy with covariance matrix adaptation. In: Proceedings of the 11th Annual conference on Genetic and evolutionary computation; 2009. p. 1411-1418.

74. Penas DR, Banga JR, González P, Doallo R. Enhanced parallel Differential Evolution algorithm for problems in computational systems biology. Applied Soft Computing. 2015; 33:86-99. https://doi.org/10.1016/j.asoc.2015.04.025