Abstract
Electroporation is a well-established technique that induces transient pores in cell membranes through intense electric fields. While previous studies have focused primarily on lipid bilayers, the role of membrane proteins in electroporation remains poorly understood. This study investigates the impact of high-intensity electric fields on the TRPV4 ion channel using molecular dynamics simulations to elucidate protein involvement in electroporation processes. The research examines two conformational states of human TRPV4 (hTRPV4, closed and inactivated) embedded in 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine lipid bilayers in a hydrated system under electric (E-)field exposure of 50, 70, and 80 MV/m. A novel algorithm was developed to precisely identify and analyze two distinct water populations: interfacial water near membrane surfaces and channel-confined water within the TRPV4 channel. Key findings reveal fundamentally different responses between these water populations. These findings provide new mechanistic insights into protein-mediated electroporation, highlighting its potential role in electroporation-driven membrane permeabilization, and suggest opportunities for developing targeted therapeutic strategies that exploit the presence and functional state of specific transmembrane proteins in membrane permeabilization processes.
Introduction
The interaction between intense electric (E-)fields and cellular membranes has been extensively investigated, particularly in electroporation, which transiently increases membrane permeability via reversible pore formation.1,2 Electroporation has demonstrated considerable potential in biomedical applications such as gene electrotransfer and electrochemotherapy, 3 driving efforts to refine models of membrane disruption, molecular transport, and intracellular dynamics. 4 In parallel, a deeper understanding of the cell pathways activated by electroporation has proven essential for optimizing its therapeutic efficacy and safety. 5
This technology is also widely used in biotechnology, 6 where it facilitates ion and molecule transport across membranes.7–9 The fundamental mechanism involves the application of high-intensity, short-duration pulsed electric field (PEF) 10 that cause transient increases in plasma membrane permeability, enabling transmembrane transport of otherwise impermeant molecules.1,10 A growing body of evidence points to the critical role of interfacial water in this process: water molecules interact with lipid headgroups, 11 helping to stabilize nascent pores and modulate the membrane’s biophysical response to E-field exposure. 12 Molecular dynamics (MD) simulations have further shown that water orientation and mobility are central to pore kinetics and electroporation mechanisms.13–16
Much of electroporation research has focused on lipid bilayers and their interplay with interfacial water; however, the cell membrane is a complex structure that also contains embedded membrane proteins, which serve critical structural and functional roles. 17 Incorporating membrane proteins into electroporation studies shifts the paradigm from lipid-only to physiologically relevant models.18,19 While PEFs are known to induce water pores, debate remains on whether ionic channels respond to nanosecond-scale electric stimuli. 20 For instance, one study 21 demonstrated that nanosecond (ns-) PEFs can activate the membrane protein phospholipid scramblase TMEM16F, indicating a potential active role for proteins in this process. Another study 20 reported the activation of the voltage-sensor domain (VSD) of the cardiac sodium channel NaV1.5 in response to nsPEF stimulation. Growing evidence indicates that voltage sensitivity extends beyond traditional ion channels to numerous transmembrane proteins in both excitable and non-excitable cells. 22 These voltage-sensitive membrane components introduce local heterogeneities that influence E-field distribution, modulating electroporation through protein–lipid interactions and localized field effects. 23 For example, in Rems et al., 24 it is shown that E-fields can induce considerably different hydration and electrostatic profiles in different voltage gated ion channels, along their transmembrane domains, resulting in the formation of complex pores stabilized by lipid headgroups. Similarly, works by this group25–27 have demonstrated that the application of intense E-fields can induce the formation of a water-filled electropore through the transmembrane domain of the aquaporin protein.
Whether protein actively drives pore formation or responds to lipid rearrangements remains an open question, but their involvement is critical for a comprehensive understanding of electroporation. This study focuses on TRPV4, a Transient Receptor Potential (TRP) ion channel broadly expressed across tissues and cell types. 28 TRP channels, particularly TRPV4, mediate osmoregulation, mechanotransduction, and thermosensation, 29 and are activated by temperature, mechanical stress, and chemical ligands.30,31 Their involvement in multiple physiological pathways and structural features makes them attractive targets for therapeutic applications 32 and suitable for investigating E-field-induced structural changes relevant to electroporation. In this study, we focused on wild-type TRPV4 channel models, reflecting the native human sequence in the absence of mutations, ligands, or artificial activators. This approach enables the receptor to be considered in physiologically relevant, ligand-free conditions, making it suitable for investigating the direct effects of E-fields on channel structure. In particular, this work investigates the role of water both at the membrane interface and within TRPV4 ion channels under the influence of intense E-fields, using an MD approach. MD simulations provide a powerful tool for analyzing conformational changes in biomolecules and their interactions with the environment, enabling detailed exploration of intense E-field exposure effects. 17 The applied field strength was carefully chosen to avoid electroporation of the lipid membrane, ensuring that observed effects primarily reflect protein-level responses.
By elucidating the interplay between interfacial water, channel-confined water, lipid dynamics, and protein behavior under high-field conditions, this research aims to advance our understanding of the fundamental mechanisms of membrane permeabilization and inform the development of electroporation-based biomedical applications.
Models and Methods
hTRPV4 ion channel modeling
To perform MD simulations, two distinct models of the TRPV4 were implemented to investigate the differential response to E-fields. Models were prepared using a standard homology modeling procedure, 33 which involves four steps: (1) template selection, (2) target-template alignment, (3) model construction, and (4) model assessment. Finally, the models have been ranked based on a quality assessment analysis. The workflow for the model implementation is illustrated in Figure 1, in the case of the closed state. In particular, the first model (“closed”) was generated using MODELLER 34 and starting from the 6BBJ.pdb structure from Xenopus tropicalis, 35 and using A0A6I8RQZ3 (Xenopus tropicalis) and Q9HBA0 (human) as template sequence and target sequence, respectively. The best model was selected from a pool of 50 output models based on the lowest Discrete Optimized Protein Energy (DOPE) scoring, which represents the most reliable parameter for identifying native-like structures among comparative models. The second implemented model (“inactivated”) followed the same workflow but was based on the 7AA5.pdb (hTRPV4). This structure reveals novel conformational rearrangements that result in an open conformation only in the presence of 4α-phorbol 12,13-didecanoate (4α-PDD), reverting to an inactivated state in its absence. The hTRPV4 “closed” and “inactivated” conformational states were of particular interest in this study, as they enclose a significantly smaller number of water molecules within the pore compared with the open state. These configurations allow us to investigate the effects of the applied E-field on a selected number of water molecules, providing clearer insights into field-induced molecular responses under low hydration conditions.

Workflow for the generation of the human TRPV4 closed state homology model using MODELLER. The process begins with the selection of the targeting sequences, the template sequences, and the TRPV4 frog model from the Protein Data Bank. A modeling protocol is then applied using a multiple sequence alignment with 86% identity, leading to the generation of 50 structural models in MODELLER. These are evaluated by Discrete Optimized Protein Energy (DOPE) scoring, and the top-ranked model is selected.
Molecular dynamics system setup and simulation parameters
The protein models have been embedded in a hydrated lipid bilayer, and the environment was created using the web-based interface CHARMM-GUI, 36 with a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) lipid bilayer dividing two ionic solutions consisting of TIP3P 37 water molecules and NaCl ions, as depicted in Figure 2 and detailed in Table 1. The dimensions of the box containing lipids and water molecules were 18 × 18 × 21 nm³, as reported in Figure 2a, allowing for the incorporation of a single protein in the center of the lipid bilayer, avoiding any problems related to periodic boundary conditions. The protein was embedded using GROMACS v. 2021.5, 38 followed by energy minimization, equilibration, and relaxation steps to ensure system stability before E-field application. The hTRPV4 receptors (Fig. 2b, c) embedded in their hydrated membrane environment, ready for the production runs, are reported in Figure 2a. Each MD simulation was performed in the constant number of particles, pressure and temperature (NPT) ensemble, with temperature and pressure kept at 310 K and 1 bar, respectively. The study was carried out considering an external E-field, oriented along the z-axis, applied normal to the lipid bilayer surface with intensities of 50, 70, and 80 MV/m, in line with similar studies in the literature. 24 Temperature coupling was performed using the velocity-rescale thermostat 39 with a relaxation time of 0.1 ps, applied separately to water and non-water groups. Pressure coupling was applied semi-isotropically using the Berendsen barostat 40 with a relaxation time of 1 ps and a compressibility of 4.5 × 10−5 bar−1 in both directions parallel and perpendicular to the bilayer. Bond lengths involving hydrogen atoms were constrained using the LINCS algorithm 41 with a fourth-order expansion. Short-range electrostatic and van der Waals interactions were cut off at 1.0 nm and treated with a potential-shift-Verlet modifier. Long-range electrostatics were computed using the Particle Mesh Ewald (PME) method, 42 with a real-space cutoff of 1.0 nm, a relative tolerance of 10−5, and a grid spacing of 0.16 nm using fourth-order B-spline interpolation. Periodic boundary conditions were applied in all three dimensions to minimize finite-size effects. Simulations were performed with a 2 fs timestep.

System’s Compounds for TRPV4 Closed and Inactivated States
human TRPV4 (hTRPV4); POPC, 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine.
Methodology for the selection of structurally relevant water molecules
A key component of this study involved developing ad hoc algorithms for the targeted selection and analysis of specific water molecules. Unlike conventional post-processing using GROMACS, VMD, or custom Python/C++ codes, this work employs a MATLAB-based approach via the GRO2MAT interface 43 to efficiently process and analyze GROMACS output data within MATLAB.
This integration enabled precise identification and tracking of two distinct water populations: (1) interfacial water near the membrane and (2) channel-confined water within TRPV4. The flowchart for interfacial water selection is reported in Supplementary Figure S1, and for each frame, phosphatidylcholine headgroup and H2O centroids are extracted, and their distances computed. The selection criterion was set, in this example, at 0.30 nm from the lipid headgroups to ensure analysis of a limited number of molecules corresponding to the immediate interfacial region of the lipid leaflet.
For the channel-confined water selection, a more sophisticated approach was required due to the complex three-dimensional geometry of the ion channel. The selection focused on water molecules located within the channel pore, specifically in the region between key structural residues that define the channel’s functional domains. According to established literature, following Nadezhdin et al., 44 we use glycine G679 and isoleucine I715 as structural reference points for the selectivity filter and activation gate, respectively, while these residues mark the narrowest regions of the pore in cryo-EM. Figure 3 provides detailed visualization of this selection strategy. Panel (a) shows the hTRPV4 closed state embedded in the phospholipid bilayer, with G679 and I715 highlighted as reference points for identifying channel-confined water within the physiologically relevant pore region. Panels (b) and (c) demonstrate the comparative visualization capabilities between VMD (upper panel) and MATLAB (lower panel) environments, highlighting the precision achieved in identifying lipid headgroups and pore residues essential for accurate water molecule selection.

The behavior of interfacial and channel-confined water molecules in hTRPV4 was analyzed in both control and E-field exposure conditions, with particular focus on the evolution of water dipole moments along the z-axis (µz) normalized by the number of water molecules during the final 15 ns preceding electroporation-induced pore formation. This normalization allows for a direct comparison of dipole alignment independent of fluctuations in water molecule count. Statistical analysis was performed using three independent simulation replicates per condition. For each replicate, dipole data were averaged over the selected time window, and group-level differences were assessed using one-way analysis of variance (ANOVA). When statistically significant effects were detected, Tukey’s post hoc test was applied to identify pairwise differences, focusing on comparisons between control and field-exposed conditions. Additionally, changes in the channel pore radius were evaluated to assess potential conformational changes associated with intense E-field exposure.
Results
The selection algorithm for specific system components was validated for accurate molecular feature identification. Figure 3 shows the VMD–MATLAB view comparison, demonstrating the reliability of the MATLAB-based approach through consistent identification of key structural elements, such as phospholipid headgroup centroids and relevant protein residues. Figure 4 shows boxplot analysis of water molecules in interfacial and channel-confined regions over 5 ns. Panels (a) and (b) show interfacial water molecules located at increasing distances (0.29–0.35 nm) from the phospholipid headgroups in the hTRPV4 closed and inactivated states, respectively. As expected, the number of water molecules increases with distance, showing a peak around 0.34–0.35 nm, consistent with the formation of the first hydration shell at the lipid–water interface. The 0.30 nm distance was selected for dipole moment analysis, capturing the strongest water–lipid interactions. Panel (c) focuses on channel-confined water molecules, displaying a representative 3D visualization. The analysis reveals that the closed conformation (green) retains a significantly larger number of water molecules within the pore region compared with the inactivated state (purple) under both control and E-field-exposed conditions, as shown in panel (d) of Figure 4. This underlines the more restricted configuration of the inactive state, which slightly hinders water retention between residues G679 and I715.

Analysis of water molecule distribution. Number of water molecules distribution at increasing distances from the lipid headgroups in the hTRPV4 closed
The global water response to intense E-field exposure and its relationship to membrane pore formation was assessed. Figure 5a shows the time evolution of the total water dipole moment along the z-axis (µz) for the hTRPV4 closed state system comparing control and E-field conditions. Under control conditions, µz remains close to zero with low-amplitude fluctuations, reflecting randomly oriented dipoles. In contrast, E-field exposure leads to a marked increase in dipole alignment along the field direction, with the most pronounced effect observed at 80 MV/m. To contextualize these results, Figure 5b illustrates the temporal progression of pore formation under 80 MV/m E-field exposure, from an intact membrane at 0 ns to visible pore opening at 18 ns, in line with the increase in dipole moment. Additional details of the membrane destabilization observed at 25 ns, which corresponds to the sharp rise in dipole alignment with an E-field intensity of 80 MV/m shown in Figure 5a, are provided in Supplementary Figure S2b. The time evolution of the total water dipole moment evaluation for the hTRPV4 inactivated state (Supplementary Fig. S3) shows a marked increase in dipole alignment along the field direction yet at the lowest E-field intensity 50 MV/m. This result confirms that the observed changes in water dipole behavior are directly associated with the electroporation process, establishing a clear mechanistic link between molecular-level water dynamics and macroscopic membrane permeabilization.

Once the electroporation threshold has been identified with values of 80 and 50 MV/m for hTRPV4 closed and inactivated state, respectively, the channel pore profile was analyzed to investigate potential configurational changes induced by intense E-field exposure, but below the threshold for bilayer electroporation onset. The pore radius along the channel’s z-axis for hTRPV4 closed state calculated using HOLE 45 is presented in Figure 6. To account for possible temporal fluctuations, the analysis was performed using three representative time points for each condition: one at the beginning of the simulation, one at the midpoint, and one just before pore formation. The resulting profiles report the mean pore radius, with standard deviations reflecting variability across the selected time points. No significant differences were observed between exposed and unexposed conditions at 50 MV/m. However, at 70 and 80 MV/m, slight deviations in the pore radius profiles reveal a misalignment likely arising from field-induced protein fluctuations. These subtle structural perturbations, though limited in magnitude, can alter the local electrostatic environment and the hydrogen-bonding network, modulating water orientation and density within the pore.

Pore radius analysis on the hTRPV4 in the closed state. The reported analysis is computed along the channel z-axis, highlighting the protein’s response to applied E-fields compared with the control condition. Mean values and standard deviations reflect measurements taken at three representative time points during each simulation. The selectivity filter and activation gate regions are highlighted in light green and orange, respectively.
A comprehensive analysis of interfacial water behavior at multiple distances from the lipid headgroups for both hTRPV4 closed and inactivated state is provided in Supplementary Figures S4 and Figure S5. Figure 7 focuses specifically on interfacial water molecules located within 0.30 nm of phospholipid headgroups, examining their dipole moment along the z-axis (µz) over 15 ns analysis windows. Panels (a) and (b) report the µz time evolution of the hTRPV4 closed and inactivated states, respectively. Notably, both channel conformations exhibit similar behavior: interfacial water molecules show no significant changes in dipole alignment across all tested conditions. The µz values remain consistently near zero in all exposure conditions, indicating that the interfacial water population retains a random orientation even under intense E-field exposure. This lack of response suggests that interfacial water molecules are strongly anchored by interactions with lipid headgroups, preventing field-induced reorientation. The consistency of this behavior across both hTRPV4 conformational states indicates that protein configuration does not influence interfacial water dynamics, highlighting the dominant role of lipid–water interactions in this region.

Analysis of interfacial water behavior under control and exposed conditions in the presence of three high-intensity E-fields. Time evolution of the dipole moment along the z-axis (μz) over 15 ns for water molecules located within 0.30 nm of the lipid headgroups, in the hTRPV4 closed and inactivated states, shown in panels
A comparable time-resolved analysis was also performed for channel-confined water, using the same 15 ns analysis window adopted for interfacial water. The resulting µz values are provided in Supplementary Figure S6 and show a clear and progressive dipole alignment under increasing E-field intensities, highlighting the marked responsiveness of channel-confined water molecules to external E-fields. To quantify this effect more rigorously, a statistical analysis was carried out for the channel-confined water population. Figure 8 presents this analysis, with panels (a) and (b) reporting the water dipole moment (µz) for hTRPV4 closed and inactivated states, respectively. Violin plots show the full distribution of µz values under control and field-exposed conditions (50, 70, and 80 MV/m), based on three independent simulation replicates per condition. Under control conditions, both channel states display dipole distributions centered near zero, indicating largely random orientations. In contrast, E-field exposure induces significant and progressive shifts toward positive µz values in both channel states. In the closed state (Fig. 8a), increasing field intensities from 50 to 80 MV/m result in progressive shifts of the distribution peak, with the 80 MV/m condition reaching mean dipole moments close to 1.4 Debye. In the inactivated state (Fig. 8b), a similar trend is observed: all field-exposed conditions lead to significant dipole alignment, but the distributions are broader and less sharply peaked compared with the closed state, particularly at higher field intensities. For both the hTRPV4 closed and inactivated conformations, statistical comparisons between each field-exposed condition and the control revealed highly significant differences (p < 0.001), reflecting a robust and consistent field-induced dipole alignment of confined water molecules. This reflects the influence of conformational state on water dynamics within the region bounded by G679 and I715 residues, which may introduce greater variability in dipole orientation despite the presence of field-induced alignment.

One-way ANOVA with post hoc tests. Panel
The clear field-dependent response of channel-confined water, contrasted with the complete lack of response in interfacial water, demonstrates the unique accessibility of channel-confined water to external E-fields. This differential behavior suggests that protein channels may serve as preferential pathways for field-induced effects, potentially contributing to localized membrane permeabilization processes.
Discussion and Conclusions
The MD study presented in this work investigates the effects of high-intensity E-fields on the hTRPV4 ion channel, with a focus on the behavior of interfacial and channel-confined water molecules. By integrating customized water selection algorithms in MATLAB with GRO2MAT tools, the study enables precise analysis of molecular responses not accessible via conventional post-processing approaches. The main observable examined in this study is the water dipole moment, as experimental and computational studies have shown that dipole alignment plays a decisive role in the onset and progression of pore formation. In particular, Bernardi et al. 25 demonstrated that water alignment within aquaporin channels contributes to the formation of transprotein electropores under E-fields. Additional work by the group of Marracino et al.13,15 proposed a dipole-driven model for electroporation, highlighting the role of water organization near lipid membranes. In this context, a central finding is the markedly different behavior between the two examined water populations. While interfacial water shows negligible dipole alignment and remains largely unaffected by E-field exposure, water molecules confined within the hTRPV4 pore exhibit significant field-induced dipole alignment. This effect is especially pronounced under 80 MV/m—the highest field strength investigated—but is also observable at lower intensities. The trend was confirmed through statistical analysis across three independent simulation replicates per condition. As shown in Figure 8, all field-exposed conditions exhibit highly significant differences compared with the control (p < 0.001) in both the closed and inactivated states, as determined by one-way ANOVA with post hoc testing. Importantly, the time window considered in these simulations excludes the onset of lipid bilayer poration. This selective polarization response suggests that channel-confined water may act as a key mediator in protein-assisted electroporation processes. The presence of aligned dipoles within the narrow channel geometry may contribute to local field enhancement, lowering the energy barrier for water permeation and possibly initiating structural destabilization of the membrane. These findings align with and are supported by previously reported experimental observations showing that E-fields can directly affect membrane protein behavior. For instance, nsPEF-induced activation has been demonstrated in TMEM16F study 21 and in the VSD of NaV1.5, 20 indicating that E-fields may modulate the structure and function of embedded transmembrane proteins. These data reinforce the rationale for including ion channels such as TRPV4 in electroporation studies and support the hypothesis that specific protein conformations can shape the local dielectric environment and serve as active participants in field-mediated permeabilization processes. Despite the application of strong E-fields, the structural integrity of hTRPV4 remains largely preserved, with only modest variation in pore radius observed at the highest field intensity. This resilience, combined with the pronounced dipole alignment of confined water, suggests that ion channels can serve as controlled permeabilization sites during electroporation. Such protein-associated pathways offer a plausible molecular mechanism for the experimentally observed enhancement of electroporation efficiency in protein-rich membrane regions. 1 The comparative analysis between hTRPV4 closed and inactivated conformations further reveals that channel geometry modulates water occupancy and response. The closed state consistently retains more water within the conduction pathway, likely due to its wider pore geometry. These conformational differences imply that the physiological state of ion channels could modulate their susceptibility to field-induced effects, enabling state-dependent control of membrane permeabilization.
While the absolute field intensities used in our simulations (50–80 MV/m) are indeed higher than those typically applied in clinical or experimental electroporation protocols (generally ranging from 0.1 to 1 MV/m), the use of such elevated fields is a well-established strategy in MD studies. This practice arises from the intrinsic timescale limitations of atomistic simulations, which require stronger fields to trigger electroporation events within nanosecond-scale trajectories. The field strengths selected in our study are consistent with a broad body of MD literature on membrane electroporation, as reported in works by Vernier et al., 46 Marracino et al., 15 Rems et al., 47 Casciola et al., 48 Delemotte et al., 49 and Fernández-Ruiz-Garate et al. 20 In support of this modeling approach, Böckmann et al. 14 showed that MD thresholds typically exceed experimental ones by a factor of 10–20, due to differences in system size, timescales, and boundary conditions (e.g., PME), which can accelerate pore formation. Moreover, since nsPEF techniques often apply megavolt-per-meter fields, our simulated intensities remain within a relevant biophysical range.
From a broader perspective, the study offers mechanistic insight into how E-fields interact with complex membrane systems at the molecular level. The findings support the concept that electric pulses can create stable, channel-like conduction pathways, 50 with water serving as an active participant in field-induced processes. This insight supports the development of targeted electroporation strategies exploiting specific transmembrane protein states.
In this study, we aimed to isolate the biophysical contributions of E-fields without confounding effects from agonist binding or allosteric modulation considering wild-type TRPV4 channel models in the absence of ligands, mutations, or artificial activators, to investigate the direct effects of E-fields under physiologically relevant, ligand-free conditions. Future work should extend this framework to other ion channels and validate predictions experimentally. Implementing a model of the hTRPV4 open state represents a promising direction for deepening the understanding of field-induced effects on functionally distinct channel conformations. Activation of TRPV4 by 4α-PDD, a well-characterized chemical agonist, promotes stable opening of the channel by binding to the intracellular S3–S4 linker region, as shown in both electrophysiological and structural studies. Modeling this ligand-induced open state would enable exploration of water and ion behavior under maximal conductance conditions and allow for comparative analysis of field sensitivity across closed, inactivated, and open conformations. This could reveal whether open-state TRPV4 channels are more susceptible to electroporation or facilitate alternative permeation pathways, offering further insight into protein-mediated field effects.
Limitations of the Study
While this study provides new mechanistic insights into protein-mediated electroporation by analyzing the dipolar response of interfacial and channel-confined water in TRPV4 ion channels, some limitations should be noted. Rather than exhaustively characterizing the channel structure under electric fields, the aim was to isolate and describe the role of confined versus interfacial water molecules as mediators of field-induced processes, offering a novel perspective on protein-mediated electroporation.
To assess the global stability of the protein, a comparative analysis of Root Mean Square Deviation (RMSD) was performed, this being a widely used metric. However, RMSD primarily captures overall conformational deviations and does not provide information on local rearrangements of helices, loops, or β-sheets, nor on transient conformational transitions that may occur under field conditions. Thus, although preliminary data on secondary structures confirmed the absence of large-scale unfolding, a more detailed structural characterization would have required further analysis such as backbone-dihedral monitoring or side-chain rotameric state evaluation, which were beyond the scope of the present work.
Future studies could integrate these approaches across a broader conformational landscape that includes TRPV4 closed, inactivated, and open states, in order to quantify domain-specific responses and assess state-dependent susceptibility to field-induced rearrangements.
Authors’ Contributions
C.P.: Writing—original draft, data curation, visualization, methodology, software, formal analysis, conceptualization, investigation; L.C.: Methodology, software, formal analysis, investigation; I.P.: Methodology, software; P.M.: Methodology, software, formal analysis; F.D.S.: Methodology, software; M.L.: Methodology, software, formal analysis, conceptualization, investigation, supervision; F.A.: Methodology, software, formal analysis, conceptualization, investigation, supervision.
