Extension of the iSoLF implicit-solvent coarse-grained model for multicomponent lipid bilayers Open Access

All-atom (AA) molecular dynamics (MD) is a powerful computational tool for investigating various biological processes in membrane environments.^1–3 To simulate moderately large biological systems for relevant time scales, a tremendous computational cost is still necessary despite the recent advancement in both hardware and software.^4–7 One popular alternative to scale up in size and time MD simulations of biological systems is using coarse-grained (CG) models. They can reduce the number of particles significantly compared to those with AA models and thereby accelerate MD simulations to reach biologically relevant time scales. Although the atomic details are lost in CG models, essential features of biomolecular structure and dynamics are still obtained in MD simulations.⁸ It should be noted that there is more than one recipe for coarse-graining the same biomolecules.^9–12 The selection of the CG model used in MD simulations, therefore, depends on the scientific question to be addressed.

CG models have been successfully applied to MD simulations of lipid membranes,^13–16 allowing the investigation of slow biological processes that would be prohibitively expensive in traditional AA MD.^17–19 The resolution of these models can widely vary from one^20,21 to a few^22–24 and up to several particles per lipid.^25,26 The coarser the mapping, the lower the computational cost. However, this comes at the expense of structural information, becoming a trade-off between computational cost and atomistic detail. Therefore, careful consideration is required in the selection of a CG model for each simulation.

The parameters in CG models can be determined following either a bottom–up or a top–down strategy. In the former, the main objective is to reproduce the structural properties of reference data, usually obtained from atomistic MD trajectories, whereas, in the latter, the aim is to reproduce experimental measurements of properties of interest. Therefore, the CG model following the bottom–up approach might lose accuracy when using it at conditions different from their parameterization. Also, it might not be easy to reproduce structural details using the CG model following the top–down strategy. For these reasons, many CG models adopt both strategies, at least partially. The electrostatics based (ELBA) polarized force field^27,28 is parameterized following the top–down approach by reproducing with accuracy different properties of lipids, such as density, surface tension, and viscosity. The Southamerican initiative for a rapid and accurate Hamiltonian (SIRAH) model²⁹ was initially developed for proteins and afterward extended to lipids.³⁰ It reproduces the physical properties of lipid bilayers, such as the area per lipid, the membrane thickness, and the order parameters. There are models such as the surface property fitting coarse-graining (SPICA),²⁶ which follow the bottom–up and top–down strategies for the bonded and non-bonded interactions, respectively. This allows the SPICA model to reproduce, for example, the phase separation in lipid bilayers at different cholesterol concentrations. One of the most popular models, Martini,²⁵ also follows top–down and bottom–up strategies for parameterizing the interactions. It has been extended to proteins,³¹ carbohydrates,³² and other molecules.³³ Martini is one of the most mature force fields, with an ecosystem built around it to facilitate the usage for studying a wide variety of systems.

There are also CG models of lipid membranes with implicit solvent approximations where the effect of water molecules is included in the interaction potentials. The main advantage of this implicit representation is the further reduction of the computational cost, allowing the simulation of even larger systems for prolonged periods. The implicit solvent variant in the Martini force field family is called Dry Martini.³⁴ Most notably, it has been used by Pezeshkian et al. to simulate an entire mitochondrial membrane.³⁵ Currently, Dry Martini is parameterized only for lipids, with the protein force field still under development.

One recent alternative is the Implicit Solvent Lipid Force Field (iSoLF).²⁴ This lipid model is designed to have a resolution similar to the Cα protein models, making it suitable to simulate lipid-protein environments in combination with the structure-based CG models.³⁶ One lipid in iSoLF consists of a linear chain containing between 4 and 6 particles, reducing computational costs further compared to Martini or Dry Martini. However, the original version of iSoLF (iSoLFv1) has several limitations: for instance, only two phospholipids (POPC and DPPC) are available together with proteins, and electrostatic interactions are not treated explicitly but are included just as effective potentials. These might cause difficulties in increasing the different types of lipid molecules available in the CG model.

In this work, we take the first step to extend and improve the iSoLF model for the new version (iSoLFv2) by addressing existing shortcomings in iSoLFv1. In the “Methods” section, we present modifications to the force field function: the explicit treatment of electrostatic interactions, the addition of polar interactions, and the introduction of new mappings to represent different lipid molecules. The “Results” section shows the parameterization results of 30 different phospholipids. We simulate the phase behavior of mixed lipid membranes composed of DOPC and DPPC at different concentrations using iSoLFv2 and compare the results with Martini v3 and Dry Martini. Finally, we discuss the status and the limitations of the iSoLFv2 model and future directions for simulating biological membranes with various lipid molecules and membrane proteins.

In both iSoLFv1 and iSoLFv2, two-tailed lipids are represented as single-tailed molecules composed of CG particles classified as hydrophilic head and hydrophobic tail beads (Fig. 1). In iSoLFv1, there are only two different head beads, H1 and H2, corresponding to the phosphatidylcholine and the phosphate bonded to the glycerol group, respectively, for POPC and DPPC. These beads are located at the center of mass of the corresponding groups they represent. The H2 particle is located between the phosphate and glycerol groups. In iSoLFv2, electrostatic interactions are explicitly treated to increase the number of available lipid types. The H2 head particle is separated into two independent CG particles, PHO and MID, representing the phosphate and the glycerol groups, respectively. Additionally, the H1 particle representing the choline group in iSoLFv1 is renamed to CHO in iSoLFv2, and four new particles are added: ETH, GLY, SER, and PHA. The first three particles correspond to the ethanolamine, glycerol, and L-serine head groups, respectively. The last one, PHA, corresponds to the unsubstituted phosphate group of phosphate lipids. In the case of the lipid tails, either two or three beads are used depending on the length of the conforming fatty acids. For 1,2-dilauroyl and 1,2-dimyristoyl lipids, two beads are used, whereas for 1,2-dipalmitoyl, 1,2-dioleoyl, 1-palmitoyl-2-oleoyl, and 1-stearoyl-2-oleoyl lipids, three beads are used. By combining tail and head beads, 30 different lipids in total can be represented with the new mapping (Fig. 1, Table S1).

We perform AA MD simulations using GENESIS^41,42 v2.0 and the CHARMM force field.⁴³ We prepare the initial configurations for 30 single-component lipid bilayers using CHARMM-GUI’s membrane builder utility.^44–48 Each bilayer is composed of 200 lipids with 100 in each leaflet. Each system contains 20 000 water molecules (100 water molecules per lipid) and is embedded in a 150 mM solution of NaCl. Next, we equilibrate the bilayers using the protocol provided by CHARMM-GUI at three different temperatures of 303.15, 313.15, and 323.15 K and a pressure of 1.0 atm. Then, we perform production runs for 1.2 µs in the NPT ensemble at the same temperature and pressure with a time step of 0.0025 ps, saving the coordinates every 10 ps. Finally, we discard the first 200 ns and use for analysis the remaining 1.0 µs, corresponding to the last 100 000 frames of the trajectories.

We perform CG MD simulations with iSoLFv2 implemented in an in-house version of GENESIS. For the parameter optimization, we prepare the initial configuration for each of the 30 single-component lipid bilayers using GENESIS_CG_tools, a set of scripts written in the Julia language.⁴⁹ Each system consists of lipids positioned in a grid of 25 × 25 nm² such that they are close to their corresponding area per lipid at 303.15 K. Table S2 summarizes, for each lipid, the number of molecules in each leaflet. We perform a short equilibration for 6.01 × 10⁵ steps in the NPγT ensemble using Langevin dynamics at three different temperatures (303.15, 313.15, and 323.15 K), a surface tension of 0.0 dyn/cm, a salt concentration of 150 mM, a time step of 10 fs, and with friction parameters for the thermostat and barostat of 0.1 ps⁻¹. During the first 1000 simulation steps, the neighbor list is updated every two steps, and afterward, every ten steps saving coordinates every 500 steps. We use the final 1.0 × 10⁵ steps of the trajectory for calculating properties (Fig. S2).

Three lipid bilayers composed of DOPC/DPPC and DOPC/DOPS mixtures at different relative concentrations (75%/25%, 50%/50%, and 25%/75%) are prepared for CG MD simulations with iSoLFv2 using GENESIS_CG_tools (Table S2). We also compare CG MD simulations of equivalent systems using GROMACS⁵⁰ v2022 with the Dry Martini³⁴ and GROMACS v2023 with the Martini 3²⁵ force fields. We build two bilayer systems of 2048 lipids (1024 per bilayer) using CHARMM-GUI’s martini builder utility.^51,52 We equilibrate the Dry Martini and Martini 3 systems using the protocol provided by CHARMM-GUI and then perform production runs using the protocol provided on the Martini webpage. We generate trajectories of 2.5 × 10⁷, 5 × 10⁷, and 1 × 10⁸ steps using a time step of 40, 20, and 10 fs and saving the coordinates every 2.5 × 10³, 5 × 10³, and 1 × 10⁴ steps for the Dry Martini, Martini 3, and iSoLFv2 systems, respectively. Finally, we use, for analysis, the last 100 frames of each trajectory.

We parameterized the intramolecular and intermolecular interactions separately, similar to the strategy used in iSoLFv1 (Fig. 2). For the intramolecular interactions, V_Bond and V_Angle, we use the multistate iterative Boltzmann Inversion method (MS-IBI).⁶⁰ This method is an extension of the original single-state IBI method and aims to improve the quality of the fitted potentials by increasing the target data averaged over states. We use AA MD simulations of 30 lipid bilayers as reference datasets, each simulated for 1 µs at three different temperatures. We map the atoms in each lipid to the corresponding CG particles. Table S3 specifies the atom names in the CHARMM topology belonging to each CG particle. After one iteration of the MS-IBI method, we obtain the distributions for each lipid's bond and angle terms. We then fit the intramolecular interaction potentials to their corresponding distributions, obtaining parameters representing the original AA structures in various temperatures. However, for DMPA, DMPE, DPPC, and DPPG, we use only the distributions at 323.15 K, and for DMPS, the distributions at 313.15 and 323.15 K, since these lipids at lower temperatures were not in the liquid-disordered phase.

The optimization is performed in four phases. First, we start with the optimization of the parameters for 1,2-dipalmitoyl-sn-glycero-3-phosphate (DLPA). In iSoLFv2, the particles used to model this lipid are PHA, MID, DL1, and DL2. However, even with four particles, exploring the entire parameter space is not computationally efficient. Thus, we reduce the dimensionality by setting ɛ_PHA = ɛ_MID and σ_PHA = σ_MID for the head beads and ɛ_DL1 = ɛ_DL2, σ_DL1 = σ_DL2, and ω_DL1 = ω_DL2 for the tail beads. Then, in the second phase, we independently optimize the σ and ɛ parameters for PC, PE, PG, and PS in the DLPC, DLPE, DLPG, and DLPS lipids while keeping all the other parameters constant and by setting the parameters of PHO equal to the ones of PHA. In the third phase, we relax the constraints for the head beads and optimize the parameters for PHO and MID independently while keeping all the other parameters constant. Finally, in the fourth phase, we keep the parameters for the head beads constant and fine-tune the parameters for the tail beads. Since all the DL lipids share the same parameters for the tails, the tuning of ɛ, σ, and ω values is guided by the minimization of the global score function, calculated as the sum of the individual score function for each DL lipid. This ensures that the final parameters for the tail beads best reproduce the properties across different lipids. For the other lipids, we use the tuned parameters for the head beads and start from the fourth phase.

We perform AA MD simulations of single-component lipid bilayers for each one of the 30 lipids considered in this work at three different temperatures. The resulting trajectories are used to parameterize both the bonded and nonbonded interactions of the CG model. In the original iSoLFv1 model, we used experimental values for fitting the nonbonded interactions. However, due to the limited experimental measurements available in the literature at different temperatures, all the reference values used in the parameterization of iSoLFv2 are obtained from atomistic MD simulation trajectories using GENESIS. Tables S4 and S5 summarize the obtained values for the APL and PPD, respectively. These values are consistent with reported experimental measurements and previous calculations with other MD software.^43,61–66 Additionally, we observe that the APL increases and the PPD decreases by increasing the temperature from 303.15 to 323.15 K, which is expected for both properties.

Bonded interactions (bond and angle) are parameterized using the MS-IBI method to match the distributions obtained from our AA MD simulations. In the iSoLF model, two-tailed lipids are represented as single-chain molecules. Thus, lipid conformations where the two tails are partially extended along the bilayer plane cannot be represented. This becomes evident when observing the Boltzmann inverted potentials for the bond interactions. Figure S3 shows these potentials for the five bond interactions in the DOPC lipid. From these interactions, only the third, fourth, and fifth bonds involve particles representing the lipid tails (DO1, DO2, or DO3). As a result, the potentials present long tails on the left side, which extend until values close to 0 Å, corresponding to structures where two consecutive beads are almost superimposed. These structures cannot be represented using a harmonic potential and, thus, represent a limitation of the model. A similar explanation follows for the Boltzmann inverted potentials for the angle interaction (Fig. S4).

For tuning the nonbonded interactions, we optimize a cost function by gradually adjusting the parameters for V_Repulsion, V_Attraction, and V_Polar potentials until the APL and PPD at different temperatures are reproduced. In our initial score function, S in Eq. (12), we only considered the L1 distances between the target and the calculated APL and PPD values. This score function guided the optimization toward parameter sets that produced stable lipid bilayers that reproduced the target APL and PPD values. However, we found that for some lipids, the resulting parameter sets favored the formation of stable pores in the membranes. These pores, being comparably small relative to the simulation box size, introduce minor deviations in the APL but have almost no effect on the PPD and cannot be detected by S. To avoid reaching those parameter sets after the optimization process, we include an additional term, H, to our score function, and name the new function S* in Eq. (13). H is defined such that it has a value of 0 when there is no pore in the membrane and a value of 1 otherwise and acts as a regularization factor during the optimization process. With this new function, we could find parameter sets that produce pore-free lipid bilayers (Fig. 3). Interestingly, for 1,2-dioleoyl lipids, the optimization process guided S* toward small values of σ_DP in the parameter space. Even though pores are not formed during the optimization, parameter sets with $σDP<7.1Å$ do not destabilize pores when used in pore-containing membranes. Therefore, only for 1,2-dioleoyl, we set $σDP=7.1Å$⁠.

Finally, using the function S*, we obtain parameter sets for all 30 lipids. An example of a complete parameterization is shown in Fig. S5. Compared with iSoLFv1, these new parameters reproduce the APL and PPD more faithfully in a broader range of temperatures (Fig. 4). It is interesting to note that lipids presenting a phase transition between 30 and 50 °C reproduce these properties partially, showing that our current model still cannot reproduce accurately transition temperatures (Figs. S6–S10). A list of all the parameters for the bonded and nonbonded interactions can be found in Tables S6–S8, with a summary of the interactions between different particles in Table S9.

To test the quality of the parameters in iSoLFv2, we performed CG MD simulations of multicomponent lipid bilayers composed of DOPC and DPPC at different concentrations. The phase behavior of these multicomponent lipid bilayers has been studied experimentally in the presence of cholesterols by Veatch and Keller.⁶⁷ In the absence of cholesterols, three zones are identified in the ternary phase diagram. At low concentrations of DPPC, it is completely dissolved in the DOPC phase, giving a single liquid-disordered phase (L_d). By increasing the concentration of DPPC at around 25%, a gel phase starts to form, producing the coexistence of both liquid-disordered and gel phases (L_d/L_β). This coexistence phase continues until DPPC reaches a threshold value of around 75%, after which a single gel phase (L_β) remains.

In CG simulations with iSoLFv2, we track the phase behavior of each component in the lipid bilayer by observing the distribution of the order parameter for each lipid type. DOPC in a single component membrane naturally occurs in a L_d phase, presenting a peak of the order parameter around 0.55. On the other hand, DPPC stays in the L_β phase and presents a peak of around 0.88. In multicomponent DOPC/DPPC simulations, we investigate the phase behaviors as a function of the relative concentrations. At a low DPPC concentration, the peaks of DOPC and DPPC are located around 0.56 and 0.66, respectively [Fig. 5(a)], indicating that both lipids are in the L_d phase [Fig. 5(d)]. In the second simulation, the membrane is in an equimolar concentration of DPPC and DOPC. We observe a shift in the peak of DPPC to the right and a spread of the distribution for DOPC [Fig. 5(b)]. This corresponds to DPPC lipids mainly in the L_β and DOPC lipids in both L_β and L_d phases [Fig. 5(e)]. Finally, in our third system, we increase the concentration of DPPC lipids, effectively inverting the concentration in our initial simulation. In this case, we observe the appearance of peaks for DOPC and DPPC around 0.86 and 0.82, respectively [Fig. 5(c)]. The location of these peaks shows that both DOPC and DPPC are in the L_β phase [Fig. 5(f)]. Figure 6 shows the three DOPC/DPPC systems' top and side views of the last frame in each trajectory.

We further investigate the behavior of lipid mixtures by comparing the neighborhood of each lipid species across different coarse-grained models. Figure 7 shows the calculated distributions of the number of neighbors for DPPC/DOPC systems at the same three relative concentrations previously described using iSoLFv2, Martini 3, and Dry Martini. These lipids share the same head group (PC) but differ in the tail particle types. In each distribution, the neighbor count is performed over lipids of the same type. Thus, when either DPPC or DOPC presented a low concentration, the peak in each distribution is shifted to the left. For the case where both lipids are present with the same concentration, the peaks are centered around the same neighbor count. A comparison in analogous systems composed of DOPS/DOPC lipids is provided in Fig. S11. In these systems, the lipid components share the same tail particle types (DO1, DO2, and DO3) but differ in the head group and present at the three concentrations a similar trend to the DPPC/DOPC systems. We also evaluate the elastic properties by comparing the bending modulus of the previously described lipid mixtures. Table I and Table S10 show the estimated values from the power spectra (Figs. S12 and S13) for the DPPC/DOPC and DOPS/DOPC systems, respectively, using the three CG models. We observe that both iSoLF and Martini 3 follow an increasing trend in the bending modulus, whereas values obtained with the Dry Martini model remain constant. Finally, we calculated the mixing entropy as described in the supplementary methods for the DPPC/DOPC and DOPS/DOPC systems, observing the phase coexistence in the DPPC/DOPC system (Fig. S14). Overall, we find that the neighborhoods are consistent among different models, and the concentration dependence of the bending modulus in Martini 3 is qualitatively reproduced in iSoLFv2.

One of the main limitations in iSoLFv1 is the implicit representation of charge interactions. This limitation is especially significant when dealing with protein–lipid interactions and when trying to capture the properties of lipid head groups with different net charges. Thus, one of the key improvements in the iSoLFv2 model is the explicit treatment of electrostatic interactions using the Debye–Hückel model, allowing the representation of charged lipids such as phosphatidyl-L-serine and phosphatidylglycerol. Additionally, the Debye–Hückel model considers the effect of the solvent by screening the electrostatic interactions as a function of the temperature of the solvent and the ionic concentration, opening the possibility to capture the dependency of bilayer properties on the salt concentration by further tuning the parameters. It is important to mention that since the dielectric constant inside lipid bilayers is comparably different from aqueous environments, we might need to refine the electrostatic interactions further, and it represents a limitation of the implicit solvent nature of the model. In the iSoLFv1 model, the center of mass of the H2 bead is localized between the phosphate and glycerol group, adding a non-zero charge to this particle would introduce artifacts due to the localization of the charge. Thus, the iSoLFv2 model represents the phosphate and glycerol with individual beads and assigns integer charges of −1e, 0e, and +1e to the negatively, neutral, and positively charged particles, respectively (Fig. 1). These integer charges are fixed and do not change during the parameterization of the non-bonded interactions. Interestingly, despite presenting charged groups in the atomic structure, the SER head group is represented with a neutral bead in the iSoLFv2 model. Thus, if we exclude the polar interaction, the SER lipid head would only interact with other particles with the repulsive potential. In fact, in the early stages of this work, we did not consider such polar interaction, making the optimization algorithm lead us to a trivial solution for the parameter set of the excluded volume interaction of the SER bead. This can be explained by analyzing the structures, APLs, and PPDs of 1,2-dipalmitoyl-sn-glycero-3-phosphate (DLPA) and 1,2-dipalmitoyl-sn-glycero-3-phospho-L-serine (DLPS). Despite differing structurally only in the presence or absence of the SER lipid head bead, both lipids present similar values for the APL and PPD and keep this trend at different temperatures (Tables S4 and S5). Since the similar property values produce similar scores using Eq. (13), the repulsive parameters for PS always converged toward zero, effectively reducing the lateral repulsion between DLPS lipids and making DLPA and DLPS the same lipid. Therefore, by adding the polar interaction between SER and charged groups, DLPS lipids could approach the target APL and PPD not necessarily by decreasing the repulsion but also by increasing the attraction between polar and charge particles, resulting in a non-trivial solution.

We find that our current parameterization best reproduces the properties of lipids that do not present a phase transition from L_β to L_d in the temperature range of 30–50 °C, as we observe by the deviations in the APL and PPD when comparing CG and AA simulations. This happens mainly when the phase transition temperature was found in the parameterization range of 30–50 °C since our model probably cannot correctly capture both phases. Additionally, the spontaneous formation of pores in membranes is addressed in iSoLFv2 by including a regularization term in the score function, as opposed to fixing the ratio between the head and tail particles as in iSoLFv1. Replacing this regularization term with explicit measurements of the line tension^68,69 will improve the score function's quality for guiding the force field's parameterization.

Finally, by reproducing the properties of lipids used routinely in experiments at different temperatures with the new parameter set, we can apply the iSoLFv2 model to a wider range of lipid environments.^70,71 However, key components are still missing in the iSoLFv2 model. The first one is the representation of sterols, specifically cholesterol. In its current state, we can represent only lipid membranes with a low to zero cholesterol content. These molecules are involved in various biological processes^72–75 and regulate lipid bilayers' mechanical and chemical properties. The second missing component is the representation of sugars. Sugars are involved in several signaling processes and including them in the model unlocks the representation of glycolipids by combining them with the current parameter set.⁷⁶ Finally, the third and most important component is the reparameterization of the lipid–protein interactions with the new parameter set since one of the key features of the iSoLF model is compatibility with other implicit solvent CG models. Decoupling the electrostatic interaction from the current, lipid–protein interaction potential is required before iSoLFv2 can be used with other protein models, and it will allow the simulation of proteins with different protonation states, something that is not possible with the current iSoLFv1 model. We are currently working on this last point, adjusting the current parameter set³⁶ to reproduce the orientation and dynamics of a wide variety of membrane proteins, and we expect soon to be able to use iSoLFv2 not only alongside the AICG2+ protein model⁷⁷ but also with other models for studying intrinsically disordered proteins⁷⁸ and liquid–liquid phase separation.^37,79,80

In this work, we have improved and extended the iSoLF implicit-solvent lipid force field, addressing various shortcomings. We have added an explicit treatment for the electrostatic interactions between charged particles and polar interactions between charged and neutral polar head particles into the force field. Additionally, we have modified the functional form for the attractive interaction between the hydrophobic tails by replacing it with a polynomial switch function, decreasing the computational cost associated with this interaction. To better represent the charge interaction between particles, we have changed the CG mapping from AA structures, representing the phosphate and the glycerol group in each lipid with individual particles instead of the one-particle representation in the original iSoLF model. By fitting the parameters to 30 different lipids, the force field can reproduce the area per lipid and the membrane thickness calculated in AA simulations at three different temperatures. Finally, we have examined the behavior of lipid bilayers composed of DOPC and DPPC at different concentrations and compared it with similar systems using the Martini 3 and Dry Martini force fields. Together, these improvements allow the iSoLF model to explore the behavior of large multicomponent systems. It will become a valuable tool once we update the interaction between lipids and proteins with the iSoLFv2 model and combine it with existing CG models. The iSoLFv2 model is publicly available in the GENESIS MD software.

See the supplementary material for additional information on the iSoLFv2 force fields, parameterizations, and comparisons between all-atom and CG MD simulations of various lipids.

This research was partially supported by RIKEN Pioneering Research Projects (Dynamic Structural Biology/Glycolipidologue Initiative) (to Y.S.), MEXT Program: Data Creation and Utilization-Type Material Research and Development Project Grant No. JPMXP1122714694 (to Y.S.), and MEXT/KAKENHI under Grant Nos. JP19H05645 and JP21H05249 (to Y.S.). We used the computer systems HOKUSAI (Project ID: Q22536, Q23536), provided by the RIKEN Information System Division, and FUGAKU (Project ID: ra00003).

The authors have no conflicts to disclose.

Diego Ugarte La Torre: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Shoji Takada: Conceptualization (supporting); Writing – review & editing (supporting). Yuji Sugita: Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.