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

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.

2 Labels and indices of sub-divided brain regions This section analyzes the estimation results in Problem 6, i.e. the stochastic and stiff VEP model at sensor level, characterized by a long seizure envelope due to a large time-scale separation.As stated in the main text, we considered a high value of noise dynamics, represented as zero-mean Gaussian noise with a standard deviation of 0.1, to achieve propagation.The result indicated that, although the system dynamics of regions corresponding to EZ and HZ are accurately estimated, the estimated trajectory for one of the PZs dampens to a stable fixed point, despite the presence of a limit cycle in the true hiddent dynamics.
In this context, we performed further tests to ascertain whether augmenting additional computational resources and refining the initial points could help us retrieve the standard parameters.We also explored a revised cost function for optimization, intending to capture more accurately the stochastic aspect of this issue.Specifically, we adopted the median of 100 simulated RMSE values for each parameter vector as the objective function.Given the computational intensity of this updated cost function, we executed these new tests utilizing SaCeSS with a notably extended time limit (58 hours) and an increased count of parallel processors (64 workers).
Table 2 displays the solutions obtained considering 4 different configurations for these tests.Each reported solution was the best out of 10 runs.The notation params true indicates the nominal value of the parameters used to generate the synthetic data.params true Same than solution 3, but starting from a different solution generated by the perturbed perturbation of params true by 10%.Figures 3 and 4 shows their results in green.
Figures 2 and 3 show the parameter values for the solutions mentioned in Table 2.In both graphs, the grey line represents the values of params true. Figure 4 illustrates the convergence of the SaCeSS runs for each test.We note that, even though the initial parameters varied significantly in quality (as indicated by cost function values), nearly all the runs converged to a common final solution in terms of the cost function.Nevertheless, as elaborated further below, these nearly identical solutions in cost do not completely align in terms of estimated parameter values.With the aim of better comparing the fitness of these different solutions, we conducted a larger set of simulations (1000) to evaluate the cost for each parameter vector, params true, solution1, solution2, solution3, and solution4.Boxplots showing the dispersion of the costs are shown in Figure 5.The median and distributions for the solutions in Table 2 are very similar, with cost function values between 2400-2500.Interestingly, this is significantly lower than the distribution provided by param true, where the median is around 3200.Therefore, although the provided solutions exhibit better fitness in terms of the RSME cost function, technically they correspond to a slight overfitting, thus explaining why we could not recover the nominal values of the parameters.However, we also suspect that there is a lack of identifiability caused by the specific cost function considered.
One well-known issue with evolutionary algorithms, metaheuristics, and other stochastic methods is that their performance can significantly vary depending on the chosen solver settings (hyperparameters).Moreover, the selection of these hyperparameters is typically problem-dependent.
In this study, we have used Optuna, a Bayesian optimization framework, to develop a methodology for automatically tuning the settings of SaCeSS in order to improve its performance on the VEP benchmarks.It is worth noting that most of the SaCeSS hyperparameters are related to the collaboration and adaptation mechanism among the cooperative workers.
We started by formulating the hyperparameter tuning as a mixed-integer optimization problem with six decision variables representing the SaCeSS configuration options to be tuned: • user options nstuck (nstuck): Maximum number of iterations for a solution in eSS without improvement before generating another solution.
• local options evalmax (evalmax): Maximum number of evaluations performed by the local solver.
• cooperation evals threshold (e thres): Parameter to determine whether to reconfigure a slave based on the number of evaluations performed without significantly improving the cooperative solution.
• cooperation mult num sendSol (mult num): Parameter to determine whether to start the reconfiguration of a slave based on the number of solutions sent by this slave.
• cooperation minimum num sendSol (min num): Minimum number of received solutions to start the reconfiguration of a slave when such slave has not sent any solution yet.
• cooperation threshold adapt (thr adapt): Value for automatic management of the reception threshold.
Table 4 details the types and ranges of these decision variables.For a more detailed explanation of these hyperparameters, please refer to the SaCeSS solver documentation.where n is the number of runs (5 in our tests), best cost i represents the fitness of the best solution reported in run i , ω is a set containing the hyperparameters of SaCeSS, which serve as the decision variables for the optimization problem, and Φ denotes static computational resource parameters related to the number of workers and the stopping time for each run i .We have selected this threshold because it is the time required to reach good quality solutions in Problem 1.
As discussed in the main text, obtaining a reasonable solution within a moderate computation time threshold can be challenging in VEP problems.For this reason, we decided to re-estimate all the problems based on the hyperparameter tuning obtained with Optuna considering Problem 1.But even for this problem this is challenging: each evaluation of the objective function involves 5 runs of SaCeSS, this operation will require at least 150 seconds, making the problem quite computationally expensive, so using standard mixed-integer nonlinear solvers would be prohibitive.
Bayesian methods are utilized in optimization when the problem size is not very large and the cost functions are expensive to evaluate, offering reasonable solutions without the need for large number of evaluations.In this context, we chose Optuna because it is widely used in similar contexts, such as in the tuning of hyperparameters of machine learning models.We employed the Bayesian optimization scheme provided by Optuna, using the default algorithm, the Tree-Structured Parzen Estimator (TPE), to determine the best SaCeSS hyperparameters for Problem 1 minimizing the objective function presented above.The Optuna solution process took 13 hours to complete, and was repeated 10 times to account for the stochastic nature of Bayesian optimization, resulting in the hyperparameters shown in Table 4.These fine-tuned settings are similar to the SaCeSS default values with the exceptions being the maximum number of evaluations for the local solver and the number of evaluations without activity before reconfiguring a worker.
April 9, 2024 7/8 Next, we applied these fine-tuned best settings to the SaCeSS solutions of the other problems.The results were very encouraging, as reported and discussed in the main text.For more details, please refer to the spreadsheet document in the supplementary information available at Zenodo https://doi.org/10.5281/zenodo.10057789.

Table 2 .
Solutions obtained by SaCeSS considering four different tests in Problem 6.The 'Initial Guesses' column indicates whether the SaCeSS optimization started from an initial parameter set.The 'Description' column provides a brief explanation of each run.solution initial guesses Description solution1 10 runs were performed using the extra computing resources mentioned random previously and optimizing the new cost function based on the median of 100 RMSE values.The parameter values of this solution are represented in red in Figure 2, with their convergence curve also shown in the same color in Figure 4. solution2 solution1 We conducted the same experiment as above, but starting from solution1 & and params true as initial points.The color blue represents the parameter params true values of this solution in Figure 2 and its convergence in Figure 4. solution3 params true 10 runs were carried out using the new cost function and the extra perturbed computational resources, starting from a solution generated by perturbing params true by 10%.Their results, in terms of parameter values and convergence, are presented in orange in Figures 3 and 4. solution4

Fig 5 .
Fig 5. Boxplot comparison of the cost function distribution for the obtained solutions.

Table 1 .
Labels and indices of sub-divided brain regions

Table 3 .
Domain of decision variables.This table presents the minimum, maximum, and default values for each variable and its type, which can be integer, float, or categorical (the latter being an integer representing a specific category).The interval column interprets the range of the variables.For example, evalmax takes various integer values multiplied by the original problem's dimension.Next, we formulated this mixed-integer optimization as the minimization a cost function defined by the geometric mean of cost values from n independent SaCeSS runs for a given problem: