AEs discover the high-variance modes instead of the slow modes

To illustrate the intrinsic limitations of classical AEs for capturing the most dominant features, we compared them and TAEs to learn a one-dimensional CV from overdamped Langevin dynamics trajectories with the following two-dimensional triple-well potential, using $ X $ and $ Y $ as the input features for NNs, setting $ \alpha $ to 1.0 and to 10.0,

(13) $$ {\displaystyle \begin{array}{l}V\left(X,Y\right)=3{\mathrm{e}}^{-{X}^2}\left({\mathrm{e}}^{-{\left(Y-\frac{1}{3}\right)}^2/\alpha }-{\mathrm{e}}^{-{\left(Y-\frac{5}{3}\right)}^2/\alpha}\right)\\ {}\hskip7em -5{\mathrm{e}}^{-{Y}^2/\alpha}\left({\mathrm{e}}^{-{\left(X-1\right)}^2}+{\mathrm{e}}^{-{\left(X+1\right)}^2}\right)\\ {}\hskip7em +0.2{X}^4+0.2{\left(Y-\frac{1}{3}\right)}^4/{\alpha}^2.\end{array}} $$

The unit of variables $ X $ and $ Y $ is angstrom, and the unit of the output of $ V\left(X,Y\right) $ is kcal mol⁻¹. The corresponding potential energy surfaces and the time evolution of the variables are shown in Fig. 2a, b, e, f. It can be easily deduced from Fig. 2a, e that for both $ \alpha = $ 1.0 and $ \alpha = $ 10.0, the important DOF should be $ X $ . As depicted in Fig. 2c, g, when training the trajectories by AEs, the value of the bottleneck layer, $ \xi $ , varies along $ X $ if $ \alpha $ is 1.0, but switches to vary along $ Y $ if $ \alpha $ increases to 10.0. This seemingly surprising result indicates that AEs erroneously learn two distinct important DOFs, when there should be only one along $ X $ . In stark contrast, TAEs behave consistently for both $ \alpha = $ 1.0 ( Fig. 2d) and $ \alpha = $ 10.0 ( Fig. 2h), and the learned variable highly correlates with the movement along $ X $ , which is also the most dominant mode. The projections of the one-dimensional learned variable onto $ X $ and $ Y $ of other models discussed in this contribution are shown in Supplementary Fig. 1. It can be rationalised by the fact that AEs actually learn the variables that correspond to the highest variance, which is, indeed, reflected in our analyses, with a variance for $ X(t) $ of 1.0, greater than that for $ Y(t) $ , equal to 0.18, when $ \alpha = $ 1.0. Conversely, when $ \alpha = $ 10.0, the variance for $ X(t) $ is 0.78, smaller than that of $ Y(t) $ , equal to 0.99. This discrepancy in the variances may also explain why AE fails to encode the correct important DOF when $ \alpha = $ 10.0. The present comparison underscores that it is crucial to introduce temporal information into the models for learning the slow modes, as opposed to the high-variance ones. Nevertheless, as discussed in a recent review of Markov state models, sometimes slow modes may still not coincide with the important DOFs that are sought (Husic and Pande, Reference Husic and Pande2018). We, therefore, need to select proper candidate CVs, or features, as the learning input, and that is why we advocate here to combine data-driven and intuition-based approaches.

Fig. 2. Potential energy surfaces of $ V(X,Y) $ with (a) α = 1.0 and (e) α = 10.0; Time evolution of X and Y when (b) α = 1.0 and (f) α = 10.0; Projections of the encoded variable $ \xi $ on X and Y from AEs training with trajectories of (c) α = 1.0 and (g) α = 10.0; Projections of the encoded variable $ \xi $ on X and Y from TAEs training with trajectories of (d) α = 1.0 and (h) α = 10.0.

Both nonlinear modified TAE and SRVs can correctly learn the slow modes

Chen et al. (Reference Chen, Sidky and Ferguson2019a) have illustrated that nonlinear TAEs learn a mixture of slow modes and high-variance modes, and proposed the nonlinear modified TAEs that can learn correctly the slow modes, but comparing to SRVs, modified TAEs cannot generate orthogonal and multidimensional latent variables that are suitable for CV-based biased simulations. Here, we confirm their results numerically by using the potential depicted in Fig. 2c, but we also propose that the modified TAE can be extended to slow feature analysis (SFA) (Wiskott and Sejnowski, Reference Wiskott and Sejnowski2002; Berkes and Wiskott, Reference Berkes and Wiskott2005; Song and Zhao, Reference Song and Zhao2022) (see Fig. 1b), which can render the latent variables orthogonal in multidimensional cases, since Eq. (4) resembles the SFA loss function in the one-dimensional case. In an SFA, if the encoder, $ \mathbf{f} $ , is multivariate, Eq. (4) should be formulated as the following generalised eigenvalue problem:

(14) $$ \Big\{\begin{array}{l}\mathbf{A}\hskip1em =\frac{1}{N-\alpha }[\hat{\mathbf{X}}(t+\tau )-\hat{\mathbf{X}}(t)]{[\hat{\mathbf{X}}(t+\tau )-\hat{\mathbf{X}}(t)]}^{\mathsf{T}}\\ {}\mathbf{A}\mathbf{V}=\mathbf{C}(t,t)\mathbf{V}\boldsymbol{\Lambda } .\end{array}\operatorname{} $$

The corresponding loss function is the sum of squared eigenvalues in Eq. (14), namely, $ \mathrm{\mathcal{L}}={\sum}_{i=1}^m{\lambda}_i^2 $ . Akin to SRVs, the CVs are then expressed as linear combinations of the components of $ \mathbf{f} $ . We note that our use of SFA in conjunction with NNs is quite similar to the deep SFA (DSFA) for change detection in images (Du et al., Reference Du, Ru, Wu and Zhang2019), except that in our two encoders, the weights and biases are shared. At first glance, Eq. (14) looks similar to Eq. (7), and also the loss function of SFA resembles that of SRVs, except that the former does not have a negative sign. Indeed, in the linear case, the two methods have been proved to be equivalent (Blaschke et al., Reference Blaschke, Berkes and Wiskott2006; Wang and Zhao, Reference Wang and Zhao2020), but the relationship is not so clear in the nonlinear case, whereby a NN with nonlinear activation functions is employed to encode $ \hat{\mathbf{X}} $ . Consequently, to investigate the performance of finding nonlinear latent variables for slow modes, in this section, we benchmark a TAE, modified TAE, SRVs and SFA using the unbiased trajectory generated from the potential energy function Eq. (13), with $ \alpha =10.0 $ . We used two hidden layers with hyperbolic tangent as the activation functions in both the encoder and the decoder part of TAEs, giving a final architecture of 2-40-40- $ n $ -40-40-2. The activation functions of the input, the bottleneck and the output layer are linear. Similarly, we used an architecture of 2-40-40- $ n $ for the other Siamese NNs (modified TAE, SFA and SRVs). We have tested both cases, where $ n=1 $ and $ n=2 $ , and used a time lag of 0.02 ps.

The projections of the encoded values in the one-dimensional case – that is, the outputs of the bottleneck layer in the TAE, the last layers in the modified TAE and the linear combinations of basis functions of SFA and SRVs – are gathered in Fig. 3a– d. In Fig. 3a, we can clearly see that the latent variable of the TAE is split into the left and right sides along $ X $ , and on each side, $ \xi $ varies along $ Y $ . This result implies that while the nonlinear TAE is able to learn a latent variable, distinguishing the two basins shown in Fig. 2c, it is also affected by the high-variance mode in $ Y $ . In stark contrast, the modified TAE, SFA and SRVs, with nonlinear activation functions, appropriately learn the slow modes, that is, the learned CV varies almost only along the transition between the two basins. In the two-dimensional case, comparing Fig. 3f, j, we note that the modified TAE learns two nearly identical modes, both varying along $ X $ . The generalisation of the modified TAE, using the SFA loss function embodied in Eq. (14), is, however, able to recognise two distinct, approximately orthogonal modes (see Fig. 3g, k), in line with the results of the SRVs (see Fig. 3h, l).

Fig. 3. Projections of the one-dimensional CV $ \xi $ learned from unbiased trajectories when $ \alpha =10.0 $ on X and Y of (a) TAE, (b) modified TAE, (c) SFA and (d) SRVs. Projections of the 2D CVs ( $ {\xi}_1 $ , $ {\xi}_2 $ ) learned from unbiased trajectories when $ \alpha =10.0 $ on X and Y of (e,i) TAE, (f,j) modified TAE, (g,k) SFA and (h,l) SRVs. Projections of the 2D CVs ( $ {\xi}_1 $ , $ {\xi}_2 $ ) learned by SRVs from an unbiased trajectory when $ \alpha =10.0 $ (m,n), from an egABF biased trajectory reweighted by Eq. (11) (o,p), and from an egABF biased trajectory reweighted by Eq. (12) (q,r).

Furthermore, we compared the two reweighting schemes, namely reweighting by uneven time intervals (Eq. (11)) and Koopman reweighting (Eq. (12)), for the calculation of the weighted correlation in SRVs from an egABF biased trajectory, with $ \alpha =1.0 $ . The variations of the two learned CVs, $ {\xi}_1 $ ( Fig. 3q) and $ {\xi}_2 $ ( Fig. 3r), using Koopman reweighting, closely correlate with the reference ones ( Fig. 3m, n) obtained from an unbiased trajectory. The results from the reweighting strategy with uneven time intervals (see Fig. 3o, p) clearly depart from the reference ones. Hence, we chose the Koopman reweighting scheme for learning CVs from biased simulations in our iterative approach. The accuracy of the scheme leaning on uneven time intervals for the evaluation of the correlations may, nevertheless, be improved by means of interpolations, either explicitly or implicitly via Fourier transform (Scargle, Reference Scargle1989), which falls beyond the scope of the present contribution.

Iterative learning of the CVs for the isomerizations of NANMA and trialanine

From the benchmarks of the previous two sections, we have concluded that SRVs are able to find both slow and orthogonal modes in multidimensional cases, and SFA performs similarly. In this section, we further test SRVs (Chen et al., Reference Chen, Sidky and Ferguson2019b) for the purpose of CV discovery based on biased simulations, applied specifically to the isomerization of NANMA and of trialanine, both in vacuum. Both peptides have been widely used as case examples in the development of novel enhanced-sampling and path-searching methods (Pan et al., Reference Pan, Sezer and Roux2008; Hénin et al., Reference Hénin, Fiorin, Chipot and Klein2010; Branduardi et al., Reference Branduardi, Bussi and Parrinello2012; Valsson and Parrinello, Reference Valsson and Parrinello2014; Tiwary and Berne, Reference Tiwary and Berne2016; Chen et al., Reference Chen, Liu, Feng, Fu, Cai, Shao and Chipot2022a). In contrast with a previous study that uses deep-TICA in a single iteration (Bonati et al., Reference Bonati, Piccini and Parrinello2021) to find from a trajectory with biasing potentials the slow mode along $ \phi $ only, we employed an iterative learning approach, akin to FEBILAE (Belkacemi et al., Reference Belkacemi, Gkeka, Lelièvre and Stoltz2022), using an initial trajectory from an egABF (Zhao et al., Reference Zhao, Fu, Lelièvre, Shao, Chipot and Cai2017) biased simulation, and performed Koopman reweighting (see Eq. (12); Wu et al., Reference Wu, Nüske, Paul, Klus, Koltai and Noé2017), as described in the Methods section. The guidelines for choosing the NN hyperparameters, the parameters of the iterative learning and the simulation details can be found in the Supplementary Material.

From the reference free-energy landscape of NANMA along backbone angles $ \phi $ and $ \psi $ ( Fig. 4a), we are able to identify the minimum free-energy pathway connecting states C_7ax and C_7eq, via the extended form, C₅ (shown as grey dots in Fig. 4b), turning to the multidimensional lowest energy (MULE) algorithm (Fu et al., Reference Fu, Chen, Wang, Chai, Shao, Cai and Chipot2020). The projection of the learned CV $ {\xi}_1 $ onto $ \phi $ and $ \psi $ in Fig. 4d clearly distinguishes the C_7ax and C_7eq states, and interprets C₅ as an intermediate state. We can also identify three basins on the PMF along the learned CV $ {\xi}_1 $ ( Fig. 4c), and by analysing the MD trajectories we have found these basins corresponding, indeed, to C_7eq, C₅ and C_7ax, respectively. Moreover, the free-energy difference between C_7eq and C_7ax obtained from Fig. 4c is equal to 2.1 kcal∙mol⁻¹, which only deviates slightly from the reference result (2.3 kcal∙mol⁻¹) deduced from Fig. 4b. The free-energy difference between C_7eq and C₅ obtained from Fig. 4c amounts to 1.2 kcal∙mol⁻¹, which is also very close to the reference result (1.0 kcal∙mol⁻¹). Additionally, the major free-energy barrier, separating C₅ from C_7ax in Fig. 4c, on the PMF along $ {\xi}_1 $ is equal to 8.2∙kcal mol⁻¹, which marginally deviates from the ground-truth reference (8.1 kcal∙mol⁻¹, marked by the red circle in Fig. 4b). In summary, not only is the learned CV, $ {\xi}_1 $ , able to discriminate qualitatively between the different metastable states encountered in the isomerization of NANMA, but the PMF along this learned variable also quantitatively predicts the correct free-energy difference and barrier.

Fig. 4. (a) Structure of NANMA and two dihedral angles $ \phi $ and $ \psi $ as candidate CVs. (b) Reference free-energy landscape of the NANMA isomerization along $ \phi $ and $ \psi $ . The grey dots show the minimum free-energy pathway (MFEP) from C_7eq to C_7ax via C₅. The dominant free-energy barrier on the MFEP is marked by the red circle. (c) The PMF along the learned CV $ {\psi}_1 $ . The three basins correspond to C_7eq, C₅ and C_7ax, respectively. (d) The value of learned CV $ {\xi}_1 $ projected on ϕ and ψ. (e) Structure of trialanine, the candidate CVs ( $ {\phi}_1 $ , $ {\psi}_1 $ , $ {\phi}_2 $ , $ {\psi}_2 $ , $ {\phi}_3 $ , $ {\psi}_3 $ ), and the eight basins (A, B, M₁, M₂, M₃, M₄, M₅, M₆) found in the free-energy landscapes along the reference CVs in (f) and the learned CVs in (g). The basins A, B, M₁, M₂, M₃, M₄, M₅ and M₆ are marked in blue, red, yellow, green, orange, white, cyan and pink, respectively. (f) Reference 3D free-energy landscape along ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ) and the corresponding MFEP (black dots). (g) 3D free-energy landscape along the learned CVs ( $ {\xi}_1,{\xi}_2,{\xi}_3 $ ) and the corresponding MFEP (black dots). (h) MFEPs found from the reference free-energy landscape along ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ) (red), the free-energy landscape along the learned CVs ( $ {\xi}_1,{\xi}_2,{\xi}_3 $ ) (blue), and the free-energy landscape reweighted from ( $ {\xi}_1,{\xi}_2,{\xi}_3 $ ) to ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ) (green).

In the paradigmatic case of NANMA, the selected candidate CVs actually coincide with the physically correct ones, namely $ \phi $ and $ \psi $ . Had we only a limited knowledge of the underlying dynamics of the process at hand, and had we included some candidate CVs that are not relevant, could our protocol still be able to learn the correct CVs? To answer this question, we tackled the more challenging example of trialanine, for which we pretend that we do not know that the three $ \phi $ dihedral angles are important (Valsson and Parrinello, Reference Valsson and Parrinello2014; Tiwary and Berne, Reference Tiwary and Berne2016) to its isomerization, and blindly select all the backbone, $ \phi $ and $ \psi $ , angles (see Fig. 4a) to form the candidate CVs, and see whether the learned CVs can render a satisfactory picture of the conformational changes. The ground-truth reference three-dimensional free-energy landscape along the three known important CVs ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ) is shown in Fig. 4f, and eight metastable states can be identified in the basins marked as A, B, M₁, M₂, M₃, M₄, M₅ and M₆. Fig. 4 depicts the three-dimensional free-energy landscape along the learned CVs ( $ {\xi}_1 $ , $ {\xi}_2 $ , $ {\xi}_3 $ ), where we can also identify eight basins. After analysis of the trajectory, we have discovered that these basins correspond to the conformations A, B and M₁–M₆, which indicates that the learned CVs are able to discriminate between the important conformations. Moreover, by applying MULE on the three-dimensional free-energy landscape along ( $ {\xi}_1 $ , $ {\xi}_2 $ , $ {\xi}_3 $ ), we determined the minimum free-energy pathway as A-M₁-M₃-B (shown as black spheres in Fig. 4g), which coincides with that found in the reference three-dimensional free-energy landscape along ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ) (shown as black spheres in Fig. 4f). It ought to be noted that a previous study (Chen et al., Reference Chen, Ogden, Pant, Cai, Tajkhorshid, Moradi, Roux and Chipot2022b) demonstrated that pathway A-M₁-M₃-B also corresponds to the most probable transition pathway (MPTP) (Pan et al., Reference Pan, Sezer and Roux2008). The free energies determined along the MFEP from the three-dimensional free-energy calculation along ( $ {\xi}_1 $ , $ {\xi}_2 $ , $ {\xi}_3 $ ) and the reference is shown in Fig. 4h in blue and red, respectively. The deviation between the blue and red curves may stem from discretization issues and difficulty to enhance sampling in the three-dimensional space (see Supplementary Material for details). However, if we reweight the free-energy landscape along ( $ {\xi}_1 $ , $ {\xi}_2 $ , $ {\xi}_3 $ ) to that along ( $ {\phi}_1 $ , $ {\phi}_2 $ , $ {\phi}_3 $ ), and identify the MFEP again, we can observe that the resulting profile (green curve in Fig. 4h) is very close to the reference one (red curve in Fig. 4h). We, therefore, conclude that the CVs ( $ {\xi}_1 $ , $ {\xi}_2 $ , $ {\xi}_3 $ ) obtained from SRVs with iterative learning reproduce the correct dynamics underlying the isomerization of trialanine, even if some fast and non-important candidate CVs are included.

Comparing the NN hyperparameters used in the NANMA and the trialanine test cases as shown in Supplementary Table 1, we can see that more neurons or computational units are used in the latter case. Therefore, if the presented approach is applied to larger biological objects, it is likely that the choice of candidate CVs differs from those for NANMA and trialanine, and the NNs will become more complex through (a) an increase of the number of layers and neurons, and (b) combination of different types of layers, for example, dropout and convolution layers, for training.

Potential connections between TAE, modified TAE and committor

In transition path theory (Weinan and Vanden-Eijnden, Reference Weinan and Vanden-Eijnden2010), considering that a molecular process characterised by two metastable states $ A $ and $ B $ , the net forward reactive flux from state $ A $ to state $ B $ can be expressed as a time-correlation function (Krivov, Reference Krivov2021; Roux, Reference Roux2021),

(15) $$ {J}_{AB}=\frac{1}{2\tau}\left\langle {\left[q\left(\mathbf{s}\left(t+\tau \right)\right)-q\left(\mathbf{s}(t)\right)\right]}^2\right\rangle, $$

where the committor, $ q\left(\mathbf{s}\right) $ , is the sum of the probability over all paths initiating from $ \mathbf{s} $ that ultimately reach $ B $ before visiting $ A $ . By definition, $ q\left({\mathbf{s}}_A\right) $ and $ q\left({\mathbf{s}}_B\right) $ are equal to 1 and 0, respectively, and as a sum of probability, $ q $ should be bounded, namely $ q\in \left[0,1\right] $ . Similar to the minimization in Eq. (4) for learning CVs, it can be envisioned that the committor function $ q\left(\mathbf{s}\right) $ can also be obtained by the minimization of Eq. (15) with the following restraints,

(16) $$ q\left(\mathbf{s}(t)\right)=\left\{\begin{array}{ll}0.0,& \mathbf{s}(t)\in A\\ {}{f}_{\theta}\left(\mathbf{s}(t)\right),& \mathbf{s}(t)\notin \left(A\cup B\right)\\ {}1.0,& \mathbf{s}(t)\in B\end{array}\right., $$

where $ {f}_{\theta}\left(\mathbf{s}(t)\right) $ is the output of a one-dimensional NN-based function with $ \theta $ as its parameters. Based on such a minimization, He et al. (Reference He, Chipot and Roux2022) have recently proposed the committor-consistent variational string (CCVS) method, where $ {f}_{\theta}\left(\mathbf{s}(t)\right) $ is parametrized by a linear combination of basis functions constructed from Voronoi cells supported by images of strings, and then optimised by an iterative Monte-Carlo procedure, as a way to determine the underlying transition pathway. At first glance, Eq. (15) looks similar to the numerator in Eq. (4) multiplied by a constant in one-dimensional cases. The main difference is that Eq. (15) is not scaled by the variance of $ {f}_{\theta}\left(\mathbf{s}(t)\right) $ or $ q\left(\mathbf{s}(t)\right) $ . Instead, as shown in Eq. (16), the committor-based loss function explicitly requires boundary conditions to identify the two metastable states, A and B, before training. In stark contrast, the loss function of modified TAEs in Eq. (4) does not feature boundary conditions, and works as a method for blind separation. The CCVS authors have further demonstrated that their method is sensitive to the diffusivities of the components in $ \mathbf{s} $ , for example, the anisotropic diffusivities along $ X $ and $ Y $ of the Berezhkovskii–Szabo potential ( Fig. 5a) resulting in different isocommittors (Weinan et al., Reference Weinan, Ren and Vanden-Eijnden2005).

Fig. 5. (a) Berezhkovskii–Szabo potential energy surface (Berezhkovskii and Szabo, Reference Berezhkovskii and Szabo2005). The latent variables projected onto (X, Y) learned by AEs, TAEs and modified TAEs in three diffusivity conditions: (b–d) $ {D}_x/{D}_y=0.1 $ , (e–g) $ {D}_x/{D}_y=1.0 $ and (h–j) $ {D}_x/{D}_y=10.0 $ . The AEs and TAEs are trained with a neural network architecture of 2-10-1-10-2 with linear activation functions used in all layers. The modified TAEs are trained with a 2-20-20-1 neural network, and the hyperbolic tangent is used as the activation functions for the two hidden layers with 20 computational units. The time lag for TAEs and modified TAEs is 10 steps.

As we can see, Eq. (15) resembles the loss function of modified TAEs embodied in Eq. (4) if $ \hat{\boldsymbol{\xi}} $ is one-dimensional. We found the similar manifestation of the anisotropic diffusivity in the CCVS (He et al., Reference He, Chipot and Roux2022) and in the learned CV intriguing, as it suggests the possibility of encoding kinetics information, such as diffusivities, in the methodology discussed in the present contributions. Here, we present a preliminary investigation of this hypothesis, wherein AEs, TAEs, and modified TAEs are compared, using the Brownian dynamics trajectories of anisotropic diffusivities sampled from the Berezhkovskii–Szabo potential (Berezhkovskii and Szabo, Reference Berezhkovskii and Szabo2005) as the training inputs. With the diffusivities along $ X $ and $ Y $ denoted $ {D}_x $ and $ {D}_y $ , respectively, we have examined three cases, namely $ {D}_x/{D}_y=0.1 $ , $ {D}_x/{D}_y=1.0 $ and $ {D}_x/{D}_y=10.0 $ . The results of the encoded variables $ \xi $ projected back onto the $ \left(X,Y\right) $ plane are gathered in Fig. 5b– j. Comparing Fig. 5b, e, h, we find that the results from the classical AE are invariant to the change of $ {D}_x/{D}_y $ . In stark contrast, in the cases of the TAEs (see Fig. 5c, f, i) and modified TAEs (see Fig. 5d, g, j), the learned CVs are affected by the different $ {D}_x/{D}_y $ ratios, which indicates that the time-series-based models are capable of reflecting the anisotropic diffusivity, and may have potential connections with the CCVS in the case of a two-state molecular process. Additionally, since AEs, TAEs and modified TAEs do not utilise restraints similar to those in Eq. (16), the learned CVs at a specific metastable state do not have a fixed value. In other words, for a committor $ q $ , we have $ q\left({\mathbf{s}}_A\right)=0 $ and $ q\left({\mathbf{s}}_B\right)=1 $ , but for a good learned CV $ \xi $ , we only know that $ \xi \left({\mathbf{s}}_A\right)\ne \xi \left({\mathbf{s}}_B\right) $ . The exact values of $ \xi \left({\mathbf{s}}_A\right) $ and $ \xi \left({\mathbf{s}}_B\right) $ are affected by the randomisation of the initial parameters and optimizers, which explains the colour flipping in Fig. 5h, contrasting with Fig. 5b, e.