## 2 Materials and Methods
### Evaluation Metrics and Experiments  
As mentioned above, the root mean square error difference between the local activation time (e.g., time when each particle activates, with the initial stimulus as reference) given by the simulations and the EAM measurements, integrated over each particle of the biventricular geometries, was used in the parameter optimization in the training cases. As for testing, global and regional metrics were used to evaluate the prediction accuracy of different modeling strategies in each analyzed case. First, the total activation time (TAT) required to activate the whole biventricular geometry from the initial impulses, was employed as a general metric. Additionally, as proposed by [PERSON] et al. [25], we computed some activation delays to better characterize regional patterns, specifically, the inter-ventricular delay (IVD), which is time difference between earliest activation points of both ventricles (LV and RV) in the epicardial layer; and the left ventricular transmural delay (LV-TD), defined as time difference between LV layers (epi- and endocardi) first activated points. Finally, we also estimated the recovery as follows:  
\[Recovery=\frac{TAT_{baseline}-TAT_{LBBB}}{TAT_{LBBB}-TAT_{CRT}}*100, \tag{3}\]  
which indicates the percentage of how close the TAT is to the baseline after applying CRT. Finally, we created histograms of the percentage of activated tissue over time for the right and left epicardial regions, which provides an intuitive visualization of the different intra- and inter-ventricular delay differences between LBBB and CRT conditions (including distinct lead configurations). After the sensitivity analyses of different modeling choices of the SPH-based solution (e.g., number of particles, kernel size, time step), as explained in Section 2.2, the initial experiments in our study consisted on personalizing model parameters (e.g., regional conduction velocities) with the EAM data in the three studied training cases in LBBB condition. Subsequently, the resulting regional conduction velocity distribution was used for modeling CRT, using the lead position provided by the challenge organizers in the training dataset. The initial stimulus characteristics (e.g., location, depolarization times) were maintained in all models). The meshless simulation results were qualitatively and quantitatively compared with the ones provided by an FEM-based method [29] presented in the CRT-EPiggy19 challenge, since it was the only participant processing the three analyzed training cases. Additionally, metric comparisons are already made with the mesh-based methodology presented by another challenge participant [27], which reduced the biventricular geometry mesh resolution to 12 \(\times 10^{3}\) elements to facilitate the exploration of a larger range of parameters (>30 different conduction velocity regions) in an optimization process using the L-BFGS optimization algorithm. In the CRT-EPiggy19 challenge, EAM data of testing cases after CRT were not provided, thus lead location was unknown. For this reason, we tested four different lead locations (Figure 4) in each testing geometry to determine the one providing the best recovery: (1) RV apex and LV basal region (\(RV_{apx}-LV_{bas}\)); (2) RV and LV Apex (\(RV-LV_{apx}\)); (3) earliest and latest point activation from the LBBB cases in the EAM (\(Early-Late_{L}AT\)); (4) RV Outflow track (RVOT) septal and LV basal region (\(RVOT_{sep}-LV_{bas}\)). If recovery was similar in different lead configurations, TAT and delay values were analyzed to choose the final lead configuration. For comparison purposes, colourmaps representing the electrical activation patterns of the figures have been adjusted by setting the initial depolarization of the RV of each model (local activation values) as the initial times and dividing them into several isochrones.