Charge carrier dynamics in conducting polymer PEDOT using ab initio molecular dynamics simulations Open Access

Conducting polymers are a unique class of organic materials characterized by conjugated single and double carbon bonds along the polymer backbone.^1–3 This conjugation allows electrons to delocalize across the backbone, creating a conducting path that can transport charges efficiently.⁴ High charge mobility in conducting polymers, combined with their low-cost and well-developed manufacturing,⁵ and versatility in chemical synthesis⁶ make conducting polymers ideal for innovative electronic applications such as organic light-emitting diodes,⁷ organic solar cells,⁸ organic field-effect transistors,^9,10 electrochemical transistors,¹¹ and organic thermoelectric generators.^12–14 

Significant progress in improving charge transport properties and carriers’ mobilities is possible with a better understanding of charge transport physics and control of the morphology of structure. In this regard, a lot of research efforts have been recently focused on understanding of charge transport in conducting polymers.^15,16 In particular, many studies provided a multi-scale description of electron transport with a realistic material morphology obtained from molecular dynamics simulations to calculate and analyze the charge mobility.^17–25 

The charge carriers in conducting polymers are polarons representing quasiparticles that are localized over several monomer units because of the strong electron-phonon coupling that induces geometrical distortion of the chain (from aromatic to quinoid).¹⁷ One of the central questions in the description of electron transport in conductive polymers is the calculation of the hopping rates polarons in the chains or between chains as.^26–31 For such calculations, there are several approaches, including the Marcus theory,^32,33 the Miller–Abrahams formula,³⁴ and the quantum nuclear tunneling model.³⁵ In these approaches, charge carriers are treated as states localized over the lattice that jump between the sites with a specific hopping rate. An alternative approach to charge carrier transport is based on the solution of the time-dependent Schrödinger equation,^36–39 which shows how charge carriers move along or between chains over time. The time-dependent approach assumes that charge carriers diffuse along polymeric chains, in contrast to the Marcus approach, where charges jump between the localized sites.

As an example of time-dependent approaches, Grozema et al.⁴⁰ studied the effect of torsional disorder on dynamics of intramolecular hole transport based on the tight-binding approximation. They focused on the electronic degrees of freedom while ignoring feedback effects from the nuclei. Also, numerical studies were conducted on polaron motion along single chains under the influence of an external electric field within a combined extended Hubbard model and the Su–Schrieffer–Heeger model.^36,41–43 Binder et al.⁴⁴ recently studied exciton dynamics in the presence of torsional fluctuations using the Multi-layer multi-Configuration Time Dependent Hartree (MCTDH) method. As a continuation of this work, Maiolo⁴⁵ performed intra-chain exciton dynamics using MCTDH, where the impact of both C–C stretch modes and ring torsions were also included. Berencei et al.⁴⁶ also studied polaron dynamics using the same generalized model (MCTDH), while Prodhan et al.^47,48 reported that charges are affected by short-range interactions, while excitons are influenced by long-range interactions like the Coulomb force between an electron and a hole. Moreover, Giannini et al.^49,50 developed an efficient fragment orbital-based surface hopping technique using mixed quantum–classical non-adiabatic molecular dynamics to study charge transport in short chains of small molecules. They then applied this approach to realistic nano-scale systems to investigate how excess charge carriers form a polaron. They demonstrated that the polaron moves through the crystal lattice by jumping between adjacent lattice sites and expanding by more than twice its size, as it causes local distortions of the lattice around it.

To date, the time-resolved dynamics of polarons in polymer chains has been carried out using either semi-empirical or model Hamiltonians, as mentioned above. However, a detailed description of the polaron dynamics relying on the first principle methods is currently missing. At the same time, first principle methods, such as density functional theory (DFT) de-facto became the methods of choice in theoretical material chemistry. It is also worth to stress that first principle DFT-based approaches can lead to results (e.g. for the electronic structure and optical transitions in conducting polymers^51–53) that are not just quantitatively but also qualitatively different as compared to the semiempirical ones. It is noteworthy that first principle methods combined with molecular dynamics simulations are widely used for study of biological electron transfer reactions as well as redox reactions.⁵⁴ 

In this paper, we develop and apply a computational approach utilizing ab initio methods, specifically the Car–Parrinello molecular dynamics (CPMD).⁵⁵ We apply this method to study polaron dynamics in conductive polymer PEDOT [poly(3,4 ethylenedioxythiophene)], which is one on the most important conducting polymers in organic- and bioelectronics.^56–58 We consider a case of intrachain polaron motion in one single chain, and interchain motion in two coupled chains with different degree of overlap. We also study the effect of temperature on polaron dynamics.

The organization of the article is as follows. In Sec. II, we introduce CPMD method for the description of polaron motion in polymer chains, and present all computational details. In Sec. III, we report molecular dynamics results for single and coupled chains. Finally, Sec. IV provides a conclusion.

Ab initio molecular dynamics (AIMD) is a computational method used to study the dynamics of molecules where the classical molecular dynamics describing the motion of atoms is combined with quantum mechanical approach to calculate the electronic structure.⁵⁹ Several AIMD methods are developed, including Ehrenfest molecular dynamics,⁶⁰ Born Oppenheimer molecular dynamics,⁶¹ and CPMD.^62,63 The differences between these methods are mainly how they combine classical and quantum equations to describe the behavior of atoms and electrons within the system. Note that applicability of AIMD relies on the Born–Oppenheimer approximation which can break down when the system involves electron dynamics in the exited states. In this case accounting for quantum nuclear effects is needed, which can be done using e.g. ab initio nonadiabatic molecular dynamics.⁶⁴ In this study we do not consider the excited states and the polaron always resides at the potential energy surface corresponding to the ground state. Thus, the Born–Oppenheimer approximation remains always valid, which justifies the utilization of AIMD.

In the present study we utilize the CPMD treating the electronic wave function as a dynamic variable that is updated at each time step based on the positions and velocities of the atoms. The electronic wavefunction is described with a set of auxiliary electronic degrees of freedom using a fictitious electronic mass introduced into the Lagrangian.⁵⁵ The lattice vibrations are described classically. Note that accounting for quantum effects in the lattice vibrations (zero-point fluctuations) can be important in some cases, e.g. for the bandlike charge transport in ultra-pure organic naphthalene crystal.⁶⁵ It remains to be seen whether zero-point fluctuations can be important for the case of charge transport in PEDOT. However, taking into account a rather amorphous character of PEDOT thin films⁵⁸ as opposed to an ideal ultra-pure crystal in Ref. 65 and much higher flexibility of PEDOT chains⁶⁶ as compared to naphthalene molecules, it is reasonable to expect that at the elevated temperature the classical description of the lattice vibration is fully justified. Note that at low temperatures the effect of the lattice vibration is practically negligible and the polaron structure in the PEDOT chain is completely determined by the electron-lattice interaction (see Sec. III A).

In this study we consider two systems, a single PEDOT chain, and two coupled PEDOT chains with different degrees of overlap, see Figs. 1(a), 2(a), and 2(f). The length of chains is chosen to be 12 monomer units.⁶⁷ The coupled chains are placed at an average π–π stacking distance of 3.5 Å by geometry optimization.⁶⁶ We use $Cl3−1$ as a counterion to computationally create a positive charge on PEDOT.⁵³ The $Cl3−1$ has a linear geometry and is placed at the distance ≈3.5 Å from the PEDOT chains. This distance exceeds the covalent bond length, such the ion is not chemisorbed on the chain. Note that the obtained results are generic in nature and do not dependent on the choice of a particular ion; our choice of $Cl3−1$ is motivated by a computational convenience where errors in the self-interaction correction of the DFT are small, such that charges at the counterion and chain are close to integer.⁶⁸ The simulation boxes for a single PEDOT chain and two PEDOT chains with full and partial overlaps are respectively 22 × 62 × 12 ų, 22 × 62 × 20 ų, and 22 × 101 × 20 ų.

A CP package of Quantum ESPRESSO⁶⁹ was used to perform CPMD. Running CPMD for our model involves several essential steps including (a) preparing the initial state, (b) displacing counterion to the opposite side of the system to facilitate electron motion, (c) reaching the electronic ground state for the system with the displaced ion, and (d) performing Car–Parrinello molecular dynamics. (Note that the considered model of the electron motion facilitated by the displacement of ions is appropriate for the case of the mixed ionic-electronic conducting polymers at low doping levels where the ion mobilities are higher than the hole ones, which gives rise to the hole-limited electrochemical doping as was recently demonstrated by Keene et al.⁷⁰)

In step (a) we start with the preparation of the initial state as shown in Fig. 1(a) for a single PEDOT chain and In Figs. 2(a) and 2(f) for two coupled PEDOT chains. In all configurations, the counterion is fixed at the top left, and the total energy is minimized through geometry optimization to reach the ground state. The electronic structure of the ground state of a neutral PEDOT chain and the chain containing a polaron is shown in Figs. S1(d) and S1(e) respectively. (Note that a detailed discussion of electron structure of PEDOT and related thiophene-based polymers at different doping levels can be found in Refs. 52, 53, and 58.) Then, in step (b) we purposefully displace the counterion, with the aim of investigating intra-chain and inter-chain charge dynamics in PEDOT chains. Specifically, for the case of a single chain, we relocate the counterion from the top-left to the top-right, as shown in Figs. 1(b) and 1(g). For the case of two coupled chains, we displace the counterion from the top-left to the bottom-right corners [Figs. 2(b), 2(g), and 2(k)]. When we move the counterion, the total energy of the system is changed, which disturbs and changes the equilibrium configuration reached in step (a). In step (c) we reach the new ground state using damped dynamics for electrons implemented in Quantum Espresso. After computing the initial electronic structure, Car–Parrinello molecular dynamics is performed in step (d) using the CPMD code.

The Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional,⁷¹ ultrasoft pseudopotentials with non-linear core corrections^72,73 and plane wave basis sets⁷⁴ are employed in all calculations. Note that range-separated hybrid functionals such as ωB97XD⁷⁵ have recently gained much attention for the prediction of properties of π-conjugated thiophene and thiophene-based oligomers⁷⁶ including PEDOT.^52,53 In contrast to the hybrid functional such as PBE and B3LYP that de-localize a polaron over the entire chain, ωB97XD predicts a polaron localization over several monomer units. However, in the presence of counterions, the hybrid functionals also predict the polaron localization due to a strong Coulomb interaction between the polaron and counterion.⁵³ Performance of the range separated functional ωB97XD and hybrid functional B3LYP (which is similar to PEB used in the present study) against experimental data for optical absorption was studied in Ref. 53. It was found that while both functionals provide a qualitatively similar description of the absorption spectra, B3LYP gives a better quantitative agreement with the experimental data. Thus, the ability of hybrid functionals to localize a polaron in the presence of counterions and a good agreement with the experimental data provide a justification for their utilization in the present calculations.

Tuning the simulation parameters plays an important role to ensure accurate and efficient CPMD simulations. In particular, for the chosen pseudopotentials, a kinetic energy cutoff of 100 and 1000 Ry for wavefunction and charge density, respectively, led to better convergency than the default ones (which are 40 and 400 Ry, respectively). For the canonical ensemble molecular dynamic simulations, time steps (δt) of 0.12 fs are used, and the longest simulation lasts up to around 720 fs (6000 steps). We use the Verlet algorithm with the fictitious electron mass 200 a.u. (the default value is 400 a.u.) during calculations of electron and atom propagation in the Car–Parrinello Lagrangian.⁷⁷ We use the Nose–Hoover thermostat⁷⁸ and set its oscillation frequency to 25 THz (the default value is 1 THz).

This section examines the polaron dynamics in single and coupled PEDOT chains for two different temperatures, 1 and 300 K. We also calculate the charge transfer rate between coupled chains and compare our findings with the Marcus theory.

To trace polaron motion, we study evolution of the charge density distribution and bond lengths alternation (BLA). The BLA is defined as a difference between the carbon–carbon bond length in the conjugated backbone of a PEDOT chain and the corresponding bond length of the neutral chain at T = 1 K, see Fig. S1. (Note that the use of the BLA as compared to the bond length is more convenient when the focus of the study is on the changes in the bond length, not the characteristics of bond length itself.⁵³)

Figures 1(a)–1(f) show polaron dynamics in a single PEDOT chain at 1 K, visualizing the BLA and charge density distribution. Figure 1(a) shows the initial state with the counterion positioned on the left. When the counterion is displaced to the right, the charge is immediately moved to the right, and becomes localized over monomers close to the new position of the counterion, as shown in Fig. 1(b). As time proceeds, charge density distribution remains unchanged, because of the strong Coulomb attraction of the localized charge to the counterion, see Figs. 1(b)–1(f).

Let us now discuss the temporal evolution of the BLA. At the beginning of the time evolution, the lattice distortion is located on the left side where the counterion was initially present, see Fig. 1(b). Following a brief period of time (20 fs), the atoms located on the right side, in the proximity to the counterion’s new position, are displaced from their initial position, as shown in Fig. 1(c). As more time passes [46 fs, Fig. 1(d)], the distortion on the left begins to disappear while the distortion on the right grows. Eventually, after about 200 fs, the bonds close to the counterion show significant deformation, whereas the atoms on the other side return to their equilibrium position, see Fig. 1(e). (Figure S2 depicts the temporal evolution of BLAs between 1 and 200 fs in shorter time steps.) After 200 fs the distortion of the lattice remains practically unchanged, see Fig. 1(f). Furthermore, the standard deviation in the BLA plot calculated over a relatively large interval (120 fs) in Fig. 1(f) is vanishingly small (not seen on the scale of the figure), which implies that the atomic vibrations in 1 K are negligible. It is worth stressing that in contrast to the case of motion of the charge that happens practically instantaneously, the BLA evolution is a much slower process because of the differences in the speeds of electron and nuclei motion.

Let us now investigate how the increase of the temperature to 300 K affects the polaron dynamics and charge localization. Figures 1(g)–1(k) show that the charge is rapidly localized on the monomers close to the re-positioned counterion and stay there during simulations. This indicates that the charge density distribution is mainly influenced by the counterion's position and increasing the temperature from 1 to 300 K does not have any effect on the charge localization. Let us not analyze the temporal evolution of the BLA. For the initial time interval (1 fs) Fig. 1(g) indicates the same behavior of the BLA as for the 1 K case. However, a close inspection reveals that at this time step the required temperature is not reached yet, see Fig. S3(a). We therefore will analyze the BLA evolution at longer times (t > 30 fs), when the temperature has reached 300 K. For these times manifestation of the lattice deformations in the time evolution of BLA are much more pronounced in a comparison to the 1 K case, see Figs. 1(h)–1(j). In particular, the gradual disappearance of the geometry distortion from one side and its subsequent appearance on the other side clearly seen at 1 K, is completely hidden by atomic vibrations. It is noteworthy that the obtained results are robust to the details of the initial distribution of the velocities of the atoms in a polymer chain, see Fig. S7.

In order to gain a better understanding of how temperature affects the behavior of the polaron and investigate the atomic vibrations, we conducted a CPMD simulations for a neutral PEDOT chain at two different temperatures, as illustrated in Fig. 3. The reason for studying the BLA in the neutral PEDOT chain is that in this case all atomic distortions stem from the temperature-induced vibrations and are not affected by the polaron. Comparing lattice distortions for a charged system and atomic vibrations for a neutral chain at 1 K reveals that the distortion caused by the polaron completely dominates the atomic vibrations. In contrast, at 300 K, the distortions induced by polaron (for a charged chain) and atomic vibrations (for a neutral chain) are of the same magnitude. In other words, at 300 K, the presence of atomic vibrations in the chain makes tracking the polaron distortion challenging because the former hides the latter. It is also worth noting that at an elevated temperature the bond lengths increase [solid blue line in Fig. 3(b)] indicating a thermal expansion of the chain. In order to hide atomic vibrations and discern a polaron formation, we averaged the BLAs over a time interval of 120 fs, as depicted in Fig. 1(k). [Note that the chosen time interval exceeds a period of atomic vibrations ≈20 fs corresponding to a typical vibration frequency of the carbon–carbon (C–C) bonds ≈1400 cm⁻¹;⁷⁹ temperature-induced vibration in a neutral PEDOT chain at T = 300 K is illustrated in Fig. S4.] This averaging provides an indication of the polaron formation on the right side of the PEDOT chain [cf. Figs. 1(f) and 1(k) showing a similar zig-zag type pattern of the BLA around bond number 40]. Finally, we note that our findings challenge a traditional picture of the polaron (based on semi-empirical or model Hamiltonians), where the effect of the atomic vibration on the polaron shape (i.e. on the pattern of the geometrical distortion of the chain) is often disregarded. In such a picture the polaron is represented as the well-defined pattern of the bond length distortion (from aromatic to quinoid structure) extended over several monomer units in the polymer chain and not affected by the temperature-induced atomic oscillations.

In this section, we examine polaron dynamics in two coupled PEDOT chains with partial and full overlaps. Figure 2 shows a time evolution of the charge density distribution and BLA for the coupled chains with a partial overlap at 1 K, and for the coupled chains with full overlap at 1 and 300 K. (The case of the coupled chains with a partial overlap at 300 K is presented in Fig. S5.) Charge density and BLA evolution at T = 1 K presented in Figs. 2(a)–2(j) shows that temporal dynamics of a polaron in two coupled chains share similar features with that one in a single chain. Namely, the charge is rapidly localized near the counterion’s new position due to the fast electron motion, and its position remains unchanged over time. Concerning the BLA evolution, the polaron distortion gradually disappears from its initial location and reappears near the new position of the counterion. These similarities hold true for both overlaps at 1 K. However, there are some differences in the polaron dynamics process between single and coupled chains.

A comparison of Figs. 1(e), 2(d), and 2(h) highlights the main difference between single and coupled chains. The difference lies in the duration of the polaron formation, i.e., how long it takes for the polaron distortion to disappear from one side and reappear on the other side of the chain when the counterion is moved there. For a single chain, as discussed above, it takes 200 fs for the polaron distortion to form near the new position of the counterion [Figs. 1(a)–1(e)]. For coupled chains with partial and full overlaps, however, the same process takes only 20 and 10 fs, respectively, as can be seen in Figs. 2(d) and 2(h). The longer time of a polaron formation on a single chain we relate to the fact that the distortions on the left and right ends of the same chain strongly interact with each other and form a long-lived standing wave extending over the entire chain. In contrast, for the case of two coupled chains such interaction is weak, and the distortions at the different chain ends do not “feel” each other and disappear (or forms) during characteristic time ∼20 fs comparable to the period of the atomic vibrations.

Another difference is related to the duration of the polaron formation for the cases of the partial and full chain overlaps, cf. Figs. 2(d) and 2(h). (Figure S6 depicts the temporal motion of BLAs in the interval between 2 and 20 fs in shorter time steps.) Specifically, the duration of polaron formation for the case of the partial overlap is twice as long as for the full overlap. The shorter time of the polaron formation for the full overlap can be attributed to the stronger interactions resulting from higher π–π overlaps.

An overall behavior of two chains at 300 K is consistent with the case of a single chain as shown in Figs. 2(k)–2(n) for the case of the full overlap. (The case of the partial overlap exhibits similar features and is presented in Fig. S5.) As in the case of the single chain, polaron dynamics is hidden behind thermal fluctuations [see Fig. 2(m)]. But, over a longer interval, polaron distortion close to the counterion can be discerned if the distortion is averaged over a time interval exceeding the period of atomic vibration, cf. Figs 1(f) and 1(k) showing a similar pattern of the BLA around bond number 40 for the lower chain.

The calculated time for the polaron formation τ, can be converted into the transfer rate K = 1/τ, giving K = 100 × 10¹² s⁻¹ and K = 50 × 10¹² s⁻¹ for the full and partial overlaps, respectively (corresponding to τ = 10 and 20 fs as discussed above). These values can be formally compared to the values of the Marcus transfer rates for the chains, see Sec. S1 in the supplementary material. For T = 300 K, K_Marcus = 8.1 × 10¹² s⁻¹ and K_Marcus = 0.5 × 10¹² s⁻¹ for the cases of the full and partial overlaps respectively. (Note that for T = 1 K the thermal hopping is practically suppressed, and the comparison to the Marcus rates is not possible as the Marcus equation is no longer applicable.⁸⁰) The difference between the hopping rates for our system and the Marcus transfer rates is related to the fact that the Marcus theory applies for the case of thermally activated hopping between the chains, whereas in our model the charge carrier hopping between the chains is induced by the displacement of the counterion between chain’s ends. We plan to report a comparison between the predictions of the CPMD simulations and the Marcus theory for the case of the thermally activated hopping in a separate study.

This study presents a new computational approach based on the ab initio CPMD to investigate temporal dynamics of charge carriers (polarons) in conducting polymer PEDOT. To study the polaron motion in a chain and between the chains we use a model where we first initiate a positively charged polaron (compensated by a negative counterion) at one end of the chain, and subsequently displace the counterion to the other end of the chain and trace polaron dynamics in the system following the evolution of the charge density distribution and bond lengths alternation. Such a model corresponds to the case of the hole-limited electrochemical doping in mixed ionic-electronic conducting polymers, where the ion mobility dominates the hole one at low hole densities.⁷⁰ We find that the charge is rapidly localized near the counterion’s new position due to the fast electron motion, and its position remains unchanged over time. Concerning the BLA evolution, for low temperature (T = 1 K) the polaron distortion gradually disappears from its initial location and reappears near the new position of the counterion. Interestingly, we find that for a single chain, it takes 200 fs for the polaron distortion to form near the new position of the counterion, whereas for the coupled chains with partial and full overlaps the same process takes only 20 and 10 fs, respectively, which is comparable to a typical period of the C–C bond vibrations. At the room temperature (T = 300 K), we find that the distortions induced by polaron, and atomic vibrations are of the same magnitude, which makes tracking the polaron distortion challenging because it is hidden behind the temperature-induced vibrations. Note that the traditional picture of polarons in thiophene-based polymers typically presented in the literature depicts polarons as a distortion from an aromatic to quinoid form localized over several monomer units (see e.g. Refs. 57 and 58). Our findings indicate that this picture can be rather oversimplified because the bond alternation caused by the presence of a charge can be completely obscured by the lattice deformations caused by the temperature-induced atomic vibration.

We believe that this study, based on the ab initio molecular dynamics, contributes to a deeper understanding of charge transport in conducting polymers and opens new avenues for further research. In particular, the AIMD approach might be instrumental in identifying the charge transport units needed in the multi-scale transport calculations.^20,23 (Charge transport units represent segments of polymer chains where charge carriers are trapped for long time, and the carrier mobility is calculated as polaron hopping between the transport units. Performing AIMD calculation and tracing polaron dynamics can be used to identify the segments where the polarons are trapped.) Using the AIMD technique one can directly calculate the diffusion coefficient (and therefore the mobilities) of polarons for intrachain and interchain motion by directly tracing the polaron dynamics (in a similar way as it has been done in e.g. Refs. 36 and 46 using the semi-empirical approaches). Note that such calculations would require utilization of the range-separated functionals⁷⁶ that localize polaron even in the absence of counterions. In the present study we investigated the dynamics of a single polaron on polymer chains. In many experimentally relevant cases polymers are, however, highly doped; for example, a doping level of PEDOT in the pristine form (i.e. as polymerized) is typically 33%.^56,57,58 The AIMD approach can be easily adapted to investigate polaron dynamics for the case of doped polymer systems. Also, the AIMD approach can be extended for the case n-doped polymers, which attract increasing attentions because of their potential applications in various fields, notably in polymer thermoelectrics.⁸¹ 

See the supplementary material for additional informative figures on the geometrical structure of a PEDOT chain, the electronic structure of the ground state of a PEDOT chain, the time evolution of BLA for a single PEDOT chain, the evolution of the system’s temperature, the temperature-induced vibration in a neutral PEDOT chain, the time evolution of BLA and charge density distribution during CPMD for coupled chains under full and partial overlap. Moreover, the calculation of charge carrier transfer rates based on the Marcus theory has been discussed in the supplementary material.

This work was supported by the European Commission through the Marie Skłodowska-Curie projects HORATES (Grant No. GA-955837). The computations were conducted on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC and HPC2N.

The authors have no conflicts to disclose.

Najmeh Zahabi: Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (equal); Project administration (supporting); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Glib Baryshnikov: Data curation (supporting); Formal analysis (supporting); Methodology (equal); Validation (supporting). Mathieu Linares: Formal analysis (supporting); Methodology (equal); Visualization (supporting). Igor Zozoulenko: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (supporting); Writing – review & editing (lead).

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