Secondary active transporters couple the transport of an ion species down its concentration gradient to the uphill transport of another substrate. Despite the importance of secondary active transport to multidrug resistance, metabolite transport, and nutrient acquisition, among other biological processes, the microscopic steps of the coupling mechanism are not well understood. Often, transport models illustrate coupling mechanisms through a limited number of “major” conformations or states, yet recent studies have indicated that at least some transporters violate these models. The small multidrug resistance transporter EmrE has been shown to couple proton influx to multidrug efflux via a mechanism that incorporates both “major” and “minor” conformational states and transitions. The resulting free exchange transport model includes multiple leak pathways and theoretically allows for both exchange and cotransport of ion and substrate. To better understand how coupled transport can be achieved in such a model, we numerically simulate a free-exchange model of transport to determine the step-by-step requirements for coupled transport. We find that only moderate biasing of rate constants for key transitions produce highly efficient net transport approaching a perfectly coupled, stoichiometric model. We show how a free-exchange model can enable complex phenotypes, including switching transport direction with changing environmental conditions or substrates. This research has broad implications for synthetic biology, as it demonstrates the utility of free-exchange transport models and the fine tuning required for perfectly coupled transport.

Introduction

Secondary active transporters use the energetically favorable flux of one ion down its electrochemical gradient to drive transport of a second substrate in the same (symport) or opposite (antiport) direction. Despite the importance of these integral membrane proteins for nutrient acquisition, metabolite transport, and toxin efflux, the mechanistic details of ion–substrate coupling are still not fully understood (Zhang et al., 2019; Yazaki et al., 2008; Amaral et al., 2014; Silverman, 2002). Frequently, transport mechanisms are depicted by models that include only the minimum number of states and transitions needed to explain transport (Oh and Boudker, 2018; Luo et al., 2018; Lolkema and Slotboom, 2019; Forrest et al., 2011; Quistgaard et al., 2016; Feng et al., 2012; Stein, 1986; Fig. 1). Through restricting transport to a single pathway, these models lead to perfect energetic coupling between substrate and ion. Occasional “slippage,” or transport through pathways apart from those prescribed by perfect coupling, is expected but assumed to be a relatively minor contributor to net flux. Although net transport stoichiometry is challenging to measure accurately, it has been measured via reversal potential for several transporters (Fitzgerald et al., 2017; Kermani et al., 2018; Nguitragool and Miller, 2006), and stoichiometric or near-stoichiometric transport was observed, as expected. Together with crystal structures (Chen et al., 2007; Iancu et al., 2013) that have captured a finite number of states for individual transporters and families, this has led to the widespread representation of such models in the literature.

However, it is becoming clear that deviations from ideal coupled transport may be common for some transporters, and that slippage pathways may even become the dominant transport pathway under certain conditions. For example, GlcPSe, a sugar transporter from Staphylococcus epidermidis, can perform sugar/proton symport, sugar uniport, and proton uniport based on pH conditions, indicating that it is possible for a single protein to engage in both coupled and uncoupled transport pathways (Bazzone et al., 2017). Some transporters engage in multiple transport regimes or stoichiometries based on substrate identity. Examples include Nramp, a metal ion transporter that performs proton-coupled symport of manganese and uniport of cadmium (Bozzi et al., 2019); or the classic sodium–iodide symporter, NIS, that performs 2:1 Na+/I+ symport but also 1:1 perchlorate symport (Dohán et al., 2007). Additionally, single point mutations of LacY result in a variety of leaky phenotypes while still performing net import of sugar molecules (Forrest et al., 2011; Guan and Kaback, 2006; Varela and Wilson, 1996). In these examples, slippage pathways significantly influence net flux. Thus, a mechanistic model that explicitly considers all possibilities for ion-coupled transport is required.

Our laboratory uses EmrE, a proton-coupled multidrug efflux pump from Escherichia coli (Purewal, 1991; Yerushalmi et al., 1995), to investigate the requirements for ion-coupled transport. Since the proton motive force is inwardly directed in E. coli, EmrE confers resistance to toxic polyaromatic cations through proton/drug antiport. EmrE’s small size, stability, and ease of purification make it a model system to study by NMR spectroscopy. This powerful technique can simultaneously provide structural, thermodynamic, and kinetic information, giving unprecedented insight into the transport cycle of an ion-coupled transporter (Robinson et al., 2017; Morrison et al., 2015, 2011; Thomas et al., 2018). Strikingly, NMR reveals that EmrE adopts states and performs transitions nominally “forbidden” for an antiporter: it can both bind drug and proton simultaneously (Robinson et al., 2017) and alternate access in every available state (Gayen et al., 2016; Morrison et al., 2015; Robinson et al., 2017). In light of this evidence, we proposed a “free-exchange” model of transport for EmrE that incorporates all experimentally observed states and transitions (Fig. 2; Robinson et al., 2017). The free-exchange model includes multiple leak pathways, yet WT EmrE demonstrates sufficient coupling to confer resistance to a wide array of toxic compounds in its native E. coli.

Here, we investigate the requirements for coupled transport using mass action kinetics to simulate the free-exchange model of transport (Horn and Jackson, 1972). Using a bottom-up approach, we first model the simplest implementation of free exchange, the 8-state model (Eq. 1; Fig. 3), before investigating the 10-state model (Eq. 2; Fig. 2) needed to describe EmrE transport (Robinson et al., 2017). We base our simulations on EmrE because it is one of the few secondary active transporters for which the rate constants are experimentally measured or estimated for all of the microscopic steps in the transport cycle. We demonstrate that only modest biasing of rates is necessary to achieve reasonably well-coupled transport, a result that has broad implications both for synthetic biology and for our understanding of the fundamental nature of active transport.

Materials and methods

Numerical simulations of transport

We model EmrE-mediated drug transport into/out of a virtual proteoliposome in a simulated liposomal-flux assay (Fig. 4). We use the default parameters in ODE15s in Matlab to numerically solve the following coupled nonlinear differential equations to determine evolution of transport over time for the generic 8-state model (Fig. 3), 
ddtEHext=k1EextHextk2EHextk13EHextDext+k14EHDext+k17EHintk18EHext ,
 
ddtEHint=k3EintHintk4EHintk15EHintDint+k16EHDintk17EHint+k18EHext ,
 
ddtEext=k1EextHext+k2EHextk9EextDext+k10EDext+k19Eintk20Eext,
 
ddtEint=k3EintHint+k4EHintk11EintDint+k12EDintk19Eint+k20Eext,
 
ddtEHDext=k5EDextHextk6EHDext+k13EHextDextk14EHDext+k23EHDintk24EHDext,
 
ddtEHDint=k7EDintHintk8EHDint+k15EHintDintk16EHDintk23EHDint+k24EHDext,
 
ddtEDext=k5EDextHext+k6EHDext+k9EextDextk10EDext+k21EDintk22EDext,
 
ddtEDint=k7EDintHint+k8EHDint+k11EintDintk12EDintk21EDint+k22EDext,
 
1VolFactExtddtDext=k9EextDext+k10EDextk13EHextDext+k14EHDext,
 
1VolFactIntddtDint=k11EintDint+k12EDintk15EHintDint+k16EHDint,
(1)
and the EmrE-based 10-state model (Fig. 2), 
ddtEHext=k1EextHextk2EHextk13EHextDext+k14EHDext+k17EHint,k18EHextk27EHextHext+k28EHHext
 
ddtEHint=k3EintHintk4EHintk15EHintDint+k16EHDintk17EHint+,k18EHextk29EHintHint+k30EHHint
 
ddtEext=k1EextHext+k2EHextk9EextDext+k10EDext+k19Eintk20Eext,
 
ddtEint=k3EintHint+k4EHintk11EintDint+k12EDintk19Eint+k20Eext,
 
ddtEHDext=k5EDextHextk6EHDext+k13EHextDextk14EHDext+k23EHDintk24EHDext,
 
ddtEHDint=k7EDintHintk8EHDint+k15EHintDintk16EHDintk23EHDint+k24EHDext,
 
ddtEDext=k5EDextHext+k6EHDext+k9EextDextk10EDext+k21EDintk22EDext,
 
ddtEDint=k7EDintHint+k8EHDint+k11EintDintk12EDintk21EDint+k22EDext,
 
ddtEHHext=k25EHHintk26EHHext+k27EHextHextk28EHHext,
 
ddtEHHint=k25EHHint+k26EHHext+k29EHintHintk30EHHint,
 
1VolFactExtddtDext=k9EextDext+k10EDextk13EHextDext+k14EHDext,
 
1VolFactIntddtDint=k11EintDint+k12EDintk15EHintDint+k16EHDint,
(2)
where the species are labeled as in Figs. 2 and 3 and the rate constants (k1k30) are defined in Table 1. VolFactInt and VolFactExt represent conversion factors between the units of the EmrE species (molecules/decimeters2) within the 2D surface of the liposome and the interior/exterior aqueous drug concentration (moles per liter). We ran sets of simulations using the initial conditions listed in Table 2 and the different combinations of rate constants as described for each experiment in Table 3. All simulations were run until the system reached steady state, as determined when the internal drug concentration changed by <1 in 1,000 over the final 5% of the run.
Parameters such as liposome size, protein:lipid ratio, and surface area of the membrane are based on our previously published liposomal-flux assays (Robinson et al., 2017). The initial conditions were also set to mimic those experimental transport assays, with the pH gradient set to higher [H]+int (lower pHint), providing an outward driving force of protons. We set initial drug concentrations inside/outside the liposome to be identical ([Drug]ext = [Drug]int) and then assessed net transport by monitoring which direction the drug moved in response to the pH gradient and other conditions in each simulation. Net transport is reported using the final value of the parameter Tr, 
Tr= [Drug]int[Drug]ext,
(3)
where [Drug]int and [Drug]ext represent the free drug concentrations inside/outside the virtual proteoliposome once steady state was reached. EmrE-bound drug is excluded from the Tr calculation. Tr starts at 1, and a value of Tr > 1 (< 1) indicates net movement of drug into (out of) the liposome during the simulation. Since the pH gradient drives outward H+ movement, Tr > 1 reflects net antiport of H+/drug, while Tr < 1 reflects net symport of H+/drug.

Constraints and simplifying assumptions

In buffered drug-monitored liposomal transport assays (Robinson et al., 2018) and in live cells (Bakker and Mangerich, 1981), the proton motive force remains relatively constant. Therefore, we make the simplifying assumption of infinite buffering, such that internal/external proton concentration remains constant throughout the simulation.

The model is also simplified by thermodynamic principles and assumptions we make based on experimental data for EmrE. In the 8-state (Fig. 3) and 10-state (Fig. 2) models, there are a total of 24 and 30 rate constants, respectively (Table 3). The thermodynamic cycle constrains drug and proton binding (pKa1 + pKd2 = pKd1 + pKaD), reducing the number of independent rate constants, where pKa = −log(Ka) and pKd = −log(Kd), with the relevant equilibria illustrated in Fig. 5.

We imposed three additional constraints based on experimental observations for EmrE: (1) rate constants for substrate-on and proton-on are diffusion limited and uniformly constant (Fig. 5 A; Morrison et al., 2015; Adam et al., 2007); (2) all on- and off-rate constants are identical for open-in and open-out conformations, which is true for EmrE due to its unique antiparallel topology and identical open-in and open-out structures (Fig. 5 B; Morrison et al., 2011); and (3) rate constants for alternating access are equal for open-in to open-out and open-out to open-in conformations (Fig. 5 C; Morrison et al., 2015, 2011; Gayen et al., 2016). It is important to note that these rate constants were all measured for EmrE in the absence of any proton gradient or transmembrane voltage, and the last two constraints may not reflect the behavior of EmrE in the presence of transmembrane gradients. Altogether, these thermodynamic and EmrE-based constraints reduce the total free parameters to 8 for the 8-state model and 10 for the 10-state model, simplifying this initial exploration of the free exchange model.

Investigating the requirements for coupled symport and antiport in the 8-state model

To test how much changing the value of the alternating-access rates alone could bias transport toward symport or antiport, we varied the alternating-access rate constants from 1 s−1 to 100 s−1 while holding the proton-off rate constants at 1,000 s−1 and drug-off rate constants at 1 s−1 (Table 3, column 1). We defined a ratio for alternating-access rate constants (Fig. 6 A) to determine the degree of bias toward antiport or symport, 
RAA=kAAantiportkAAsymport,
(4)
where 
kAAantiport=kAAEH=kAAED,
(5)
and 
kAAsymport=kAAE=kAAEHD.
(6)
RAA < 1 favors flux through the E and EHD states and thus symport, while RAA > 1 favors flux through the EH and ED states and thus antiport. The different states of EmrE are as defined in Figs. 2 and 3 and reflect EmrE (E) bound to proton (H) and/or drug (D). To reduce the number of parameters, a single antiport alternating access rate was used for both antiport states, and a single symport alternating access rate was used for both symport states. This coupling of the antiport and symport rates results in the best-case scenario for drug and proton coupling at a given RAA.
To test how much changing the value of the substrate-off rate constants alone could favor symport or antiport phenotypes, we varied proton-off rate constants from 1 to 1,000,000 s−1 and drug-off rate constants from 0.01 to 10,000 s−1 (such that pKa values ranged from 4.0 to 10.0 and Kd values ranged from nM to mM, respectively) while holding alternating-access rate constants uniformly constant at 1, 10, or 100 s−1 such that RAA = 1 (Table 3, column 2). To maximize coupling efficiency and reduce the dimensionality of the parameter space within the constraints of the thermodynamic cycles, we imposed koffED/koff DEHD = koffEH/koff HEHD (Fig. 6 B) and defined this ratio as Roff: 
Roff=koffEDkoff DEHD=koffEHkoff HEHD.
(7)
Roff > 1 favors flux through E and EHD states and thus symport, while Roff < 1 favors flux through EH and ED states and thus antiport.

Finally, to test whether off-rate constants and alternating-access rate constants were additive or synergistic, we varied both Roff and RAA (Table 3, column 3). Alternating-access rate constants were set at 1, 2, 3, 4, or 5 s−1 to span 0.125 ≤ RAA ≤ 5.

Modeling the effect of a second protonation event

The liposomal assay simulations were repeated using the 10-state model, which includes a second protonation event. Since protonation is sequential, this model includes an additional constraint that the first proton binds with higher affinity than the second proton, pKa1pKa2. This negative linkage between two proton binding events is expected for two protonation sites in close proximity within a hydrophobic environment, as observed for EmrE (Morrison et al., 2015). We again varied substrate-off rate constants, with simulations performed for six values of pKa2 (Table 3, column 4). All alternating-access rate constants were held constant at 1 s−1.

Investigating influence of environmental pH on transport

Our initial simulations all used an inverted-physiological pH parameter with lower pH on the inside of the membrane to mimic our liposomal flux assays (Robinson et al., 2017). However, to explore how pH affects transport under conditions more closely aligned to the environment of a transporter in the bacterial inner membrane, we performed additional simulations holding pHint = 7.4 while varying pHext ±2 units (pHext = 5.4–9.4). This was tested in both the 8-state model and 10-state model to see if the extent of drug transport drug differed depending on the number of proton binding events. For the 8-state model, we used an antiport configuration with RAA = 10 and Roff = 0.1 (Table 3, column 5); for the 10-state model, we used rates estimated for EmrE interacting with the substrate tetraphenyl phosphonium+ (Table 3, column 6; Robinson et al., 2017).

Investigating how a single transporter can act as both an antiporter and a symporter

For these simulations, we varied the drug-dependent rate constants while holding constant the protein-specific rate constants. Rate constants that vary with the identity of the transported drug include drug-off (koffED, koff DEHD), proton-off for the drug-bound state (koff HEHD), and alternating-access rate constants in the ED and EHD states (kAAED, kAAEHD). Protein-specific variables are invariant to drug identity and include alternating-access rate constant in the apo, singly, and doubly protonated states (kAAE, kAAEH, kAAEH) as well as proton-off rate constants (koffEH and koffEHH, which dictate pKa1 and pKa2, respectively). We ran simulations scanning through physiological ranges for koffED, koffEHD, kAAED, and kAAEHD for two different pKa combinations: pKa1 = 8.2 and pKa2 = 7.0, representing WT EmrE (Table 3, column 7), and pKa1 = 7.0 and pKa2 = 5.0, representing a mutant with pKa values downshifted toward more typical glutamate pKa values (Table 3, column 8). Tr was used to assess whether the downshifted mutant permitted both antiport and symport upon changes to drug-dependent rates. These experiments followed pH conditions for liposomal transport assays.

Results

Building mass-action kinetic models of free exchange

Our 8-state model includes all states required to describe both coupled antiport and coupled symport of one proton and one drug (Fig. 3, B and C), while our 10-state model includes the experimentally observed states required to describe two-proton–coupled transport by EmrE (Fig. 2). Both models contain three types of rate constants that determine net transport phenotype: substrate-on rate constants, substrate-off rate constants, and alternating-access rate constants. These are constrained three different ways based on observations of EmrE. First, we assume that (1) substrate on-rates are diffusion-limited (Adam et al., 2007; Fig. 5 A). For EmrE, all drug substrates studied in detail have similar on-rates, and the off-rate determines each binding affinity, which ranges over five orders of magnitude (Adam et al., 2007; Morrison and Henzler-Wildman, 2014). Similarly, we assume that proton on-rate is uniformly constant and diffusion limited, and off-rate determines pKa (Morrison et al., 2015; Adam et al., 2007). Next, since EmrE has an antiparallel topology that results in identical open-in and open-out conformations (Morrison et al., 2015), we assume that (2) binding affinities (and thus on- and off-rates) are identical for open-in and open-out states (Fig. 5 B) and that (3) alternating-access rates are equal in both directions (Fig. 5 C; Morrison et al., 2015, 2011; Gayen et al., 2016). As a result, substrate-off rate constants and alternating-access rate constants are the two groups of rate constants that are variable parameters in our simulations.

We model EmrE reconstituted in a virtual proteoliposome with infinite buffering both inside and outside, resulting in constant pH throughout the simulation. Each simulation is run using a particular set of rate constants (see Table 1) for the microscopic steps in the transport mechanism. The same initial conditions are used unless otherwise noted, with pH lower inside than outside to create a pH gradient that will drive protons out of the liposome (see Table 2). Each simulation continues until steady state is reached, as assessed by monitoring the change in [Drug]int over time. Different experiments comprise sets of simulations (Table 3) that probe the influence of pH, alternating-access rates, and substrate off-rates on the final drug gradient achieved at steady state. The direction and extent of the drug gradient is assessed using the ratio of Tr (Eq. 3), where Tr < 1 represents movement of drug out of the liposome (net symport) while Tr > 1 represents movement of drug into the liposome (net antiport). For perfectly coupled stoichiometric transport, the steady-state Tr is empirically calculated from the initial proton gradient and the integer coupling stoichiometry constant n (Nguitragool and Miller, 2006): 
Trmax= ([Hint]+[Hext]+)n.
(8)
In our model, Tr at steady state is determined by ΔpH and the relative values of the rate constants. It is independent of the exact value of the initial drug concentration, which affects the rate at which steady state is achieved but not the steady-state Tr. Comparing the Tr for each simulation to Trmax provides insight on how tightly coupled the net transport process is. We defined transport as highly coupled if the Tr for that simulation was at least 80% of Trmax. For example, for the 8-state model with n = 1, Trmax would be 10.0 (0.1) for antiport (symport), with highly coupled transport defined as Tr > 8 (antiport) or < 0.125 (symport).

Requirements for highly coupled transport of one H+ and one substrate

Fig. 3 depicts the antiport and symport pathways within the 8-state model. While multiple pathways exist for both modes of transport, each contain key transitions that couple proton and drug flux across the membrane: alternating access in the EH and ED states for antiport or the E and EHD states for symport. For a complete cycle of antiport, EH must alternate access to translocate bound proton in one direction, while ED must alternate access to translocate drug in the opposite direction. Similarly, a complete cycle of symport requires alternating access of EHD to simultaneously translocate bound drug and proton in one direction before alternating access of E returns the empty transporter to the opposite side. These are exactly the states where alternating access occurs in the single-pathway, tightly coupled transport models shown in Fig. 1.

If all rate constants in the free exchange model are equal, all pathways are equally likely, and the transporter simply facilitates downhill diffusion of both proton and drug. To achieve coupled transport, rate constants must be skewed to tip the balance between transport regimes by increasing the relative flux through the key transitions for antiport or symport. We hypothesized that this could be accomplished through manipulating the rates for substrate-off and/or alternating access. Varying the alternating-access rate constants will directly impact the relative frequency of these transitions that are critical for moving substrate across the membrane. On the other hand, varying substrate off-rate constants alters the partitioning between futile substrate binding and release and productive movement of molecules across the membrane via the key alternating-access transitions. We explored the independent and combined effects of substrate off-rates and alternating-access rate through three experiments to explore the transport behavior of the free-exchange model.

In the first set of simulations, alternating-access rate constants were varied around a physiological range while holding proton off-rate constants at 1,000 s−1 and drug off-rates constants at 1 s−1 for all states of the transporter (Table 3, column 1). Alternating-access rate constants were grouped according their contribution to antiport or symport and covaried to reduce the number of free variables within the model. We examined how shifting the ratio RAA (Eq. 4) of the key alternating-access rate constants (Eqs. 5 and 6) affected Tr. For values of RAA > 1, the key antiport states alternate access faster, while for RAA < 1, this is true for the key symport states. Fig. 7 A shows the relationship between RAA and Tr. To reach the 80% Trmax benchmark for highly coupled transport, the alternating-access rates need to be skewed 25-fold in either direction: RAA ≥ 25 for antiport or RAA ≤ 1/25 for symport. This reveals the expected symmetry within the 8-state free-exchange model: whatever degree of biasing is required for highly coupled antiport is the inverse of what is required for highly coupled symport.

In the second set of simulations, we varied substrate off-rate constants across the range of values observed for EmrE substrates while holding alternating-access rates constants uniformly at 1, 10, or 100 s−1, resulting in a single kAA value for all alternating access transitions (Table 3, column 2). To minimize the number of free variables needed for simulation, we varied the substrate off-rate constants critical to either antiport (koffEH and koffED) or symport (koff HEHD and koff DEHD) such that 
koffEDkoff DEHD=koffEHkoff HEHD
(Fig. 6 B). We then examined how shifting the ratio of these rate constants (Roff, Eq. 7) affected transport. For Roff < 1, the antiport-critical states have slower off-rates and higher stability than the symport-critical states, indicating competitive binding between drug and proton as expected for antiport. For Roff > 1, this is inverted, and drug and proton bind cooperatively. Similar to the RAA experiment, Tr was calculated for each value of Roff (Fig. 7 B). Although RAA = 1.0 and neither symport nor antiport is directly favored by the alternating access steps in the transport cycle, alternating-access rates do affect transport behavior when they are the same order of magnitude as drug and proton off-rates. If alternating access is fast, then alternating access is less likely to result in coupled movement of molecules across the membrane. This is particularly important for symport, where both substrates must bind on one side of the membrane and be released on the other side. Thus, at faster alternating access rates (Fig. 7 B, orange and yellow lines), the antiport and symport phenotypes are no longer symmetric (inverting Roff does not result in inverted Tr), as symport is less well coupled. At slower alternating access rates, this symmetry is restored. When kAA = 1 s−1, 80% Trmax is achieved by 630-fold off-rates bias in either direction (Roff ≤ 1/630 for antiport or Roff ≥ 630 for symport).

In the third set of simulations, we varied both Roff and RAA. Fig. 7 D shows Tr as a function of Roff similar to the second set of simulations; however, each curve now represents a different RAA value (Table 3, column 3). With fivefold skewing of alternating-access rate constants (Fig. 7 D, green line), only a 23-fold difference in the substrate off-rate constants is needed to achieve 80% Trmax, much less than the 630-fold-difference required if alternating-access rates are uniformly held at 1 s−1 (Fig. 7 B). This result demonstrates that the combined effects of alternating-access and substrate off-rates synergize nonlinearly such that robust proton-coupled transport phenotypes can be achieved without the need for highly skewed rates. In other words, highly coupled transport does not require states or pathways to be excluded from the transport model. Only relatively modest biasing of rate constants channels the majority of flux through key pathways and achieves relatively well-coupled transport.

Modeling the transport cycle of EmrE

EmrE requires a more complex transport mechanism due to its ability to perform two-proton coupled antiport (Fig. 2). Since there is no evidence that EmrE can simultaneously transport drug and two protons, the 10-state model introduces an asymmetry. Whereas the 8-state model possesses a maximum of 1 coupled proton per drug (n = 1, Eq. 8) for both antiport and symport, the 10-state model increases n for antiport to n = 2, while symport remains at n = 1. In other words, an initial 1.0 pH unit gradient can now drive a maximum 100-fold drug gradient (Tr = 100) via antiport but can still only drive a maximum 10-fold drug (Tr = 1/10) gradient via symport. In effect, this extra protonation event creates a new key antiport state, EHH (EmrE bound to two protons, see Fig. 2), that cannot contribute to symport. We explored how this asymmetry affects the distribution of transport outcomes.

Since the second protonation event introduces an additional proton off-rate (koffEHH), we repeated the experiment varying Roff (Fig. 7 B, blue trace) but over six different values of this newly introduced proton off-rate (Fig. 8). Proton on-rates are assumed to be diffusion limited and constant at 1010 s−1, so this is recapitulated as six curves of different pKa2 values for the second protonation event (Table 3, column 4). Protonation appears to be sequential in EmrE, thus requiring pKa1pKa2 (Morrison et al., 2015). This restricts the parameter space, resulting in the gray portion of Fig. 8 (for example, if pKa2 = 9, pKa1 can range only from 9 to 10). At low pKa2 values, the second protonation event is rare, and the 10-state model greatly resembles the 8-state model (Fig. 8, overlapping orange and blue lines). However, as the value of pKa2 rises above 7, the EHH state becomes populated, and transitions through this state contribute to net flux, resulting in increasing domination of antiport phenotypes. In fact, as pKa2 approaches 9, the system starts to resemble a stoichiometric transport model, where perfect competition yields perfectly coupled antiport. Thus, the second protonation event biases the transporter toward antiport under these particular conditions with only ΔpH and no membrane potential.

Effect of pH on steady-state transport phenotypes

Next, we examined the influence of environmental conditions on net transport. In our simulations, the energy required to drive coupled transport is stored in the proton gradient. For perfectly coupled stoichiometric transport, a given ΔpH will drive a set substrate gradient according to Eq. 8. However, in free-exchange transport, the degree of proton/drug coupling is inherently influenced by the exact internal and external pH values. As proton on-rates depend on both the proton on-rate constant and the proton concentration, a change in pH affects the distribution of the key states and the partitioning between different pathways. Thus, even for a constant ΔpH, a change in mean pH alters the balance of antiport, symport, and uncoupled uniport.

To investigate this phenomenon, we ran simulations in the 8-state and 10-state models, holding the rate constants while varying pH. The rate constants were chosen to favor antiport using EmrE-like rates for both models (Table 3, column 5 and 6, for the 8- and 10-state model, respectively). Internal pH was held at physiological levels (pHint = 7.4) while external pH was varied 2 units higher or lower. Fig. 9 shows how the two models diverge in behavior. Once again, the 8-state model displays symmetric behavior, with an inverted pH gradient leading to an inverted Tr. Interestingly, this is not the case for the 10-state model, which drives a much larger drug gradient with an acidic external environment (drug efflux, lefthand side of Fig. 9) than with a basic external environment (drug uptake, righthand side of Fig. 9). This reveals that a 10-state transporter is a better antiporter under physiological conditions than it is with an unnatural inverted pH gradient, providing a possible failsafe to resist backflow of drug into a cell when the cell is exposed to a basic external environment.

How an antiporter can become a symporter

The W63G point mutation of EmrE abolishes classic WT substrate specificity and replaces it with resistance to macrolide antibiotics such as erythromycin. However, it also confers concentrative uptake of aliphatic polyamines into E. coli (Brill et al., 2012). This combination of erythromycin export and polyamine import means that W63G-EmrE must switch transport direction based on substrate identity alone. While it is not experimentally resolved whether the polyamine import phenotype is the result of proton-coupled symport or of uniport with positively charged polyamines driven into the cell by the negative-inside membrane potential, here we explore how a single transporter may switch between coupled antiport and coupled symport using the 10-state model.

For EmrE, the identity of the bound drug can change the alternating-access rate over nearly two orders of magnitude (Morrison and Henzler-Wildman, 2014). Additionally, EmrE can bind drug with affinities ranging over five orders of magnitude, with drug-binding affinity determined by off-rate for those substrates that have been studied in detail (Adam et al., 2007), and the drug off-rate varies with protonation state of the transporter (Robinson et al., 2017). Thus, we completed simulations varying all drug-dependent rate constants (kAAED, kAAEHD, koffED, koff HEHD, and koff DEHD) over reasonable physiological ranges around the known values for EmrE. Since pKa2 must be <7 to allow for symport in the 10-state model, we ran two simulations (Table 3, columns 7 and 8): one with values estimated for WT EmrE (Fig. 10 A, pKa1 = 8.2 and pKa2 = 7) and one with lower pKa values (Fig. 10 B, pKa1 = 7.0, pKa2 = 5.0), closer to what is expected for a glutamate residue in aqueous solution. These figures are graphed in 2D space by collapsing alternating-access rates and drug/proton off-rates as ratios. As expected, no combination of drug-induced rates can drive symport when using values estimated for WT EmrE (Fig. 10 A) under these conditions. However, with sufficient lowering of both pKa values, coupled antiport and symport can coexist in a single transporter (Fig. 10 B). This simple scheme allows for the identity of transported drug to determine the type of proton-coupled transport, a striking result.

Discussion

Transporters have been traditionally classified as antiporters, symporters, or uniporters. Mechanistic models of transporters generally reflect these strict classifications and include only pathways that allow for stoichiometric coupled transport (Lolkema and Slotboom, 2019; Oh and Boudker, 2018; Stein, 1986). Recently, improved experimental methods applied to a broader set of transporters have revealed that some transporters cross these boundaries and do not fit cleanly within a single class (Bazzone et al., 2017; Bozzi et al., 2019; Dohán et al., 2007; Nguitragool and Miller, 2006; Robinson et al., 2017). In light of this, a new generalized model of transport is required that can simultaneously accommodate all transport modes. We initially developed an unrestricted free-exchange model to account for all the observed states and transitions of EmrE (Robinson et al., 2017). The free-exchange model theoretically allows this proton-coupled drug efflux pump to perform antiport, symport, and uniport. Here, we have shown that such an unrestricted model can describe the antiport phenotype seen in EmrE (Fig. 10 A), and in fact the second proton binding event favors antiport (Fig. 8) and enhances the efficiency of antiport relative to symport under these simple ΔpH-only conditions (Fig. 9). These simulations can also recapitulate the curious phenotype of W63G (Brill et al., 2012), a point mutant of EmrE in which transport direction reverses based on the identity of the transported drug (Fig. 10 B).

Most importantly, we have shown that highly coupled transport can be achieved in the free-exchange model with minimal biasing of rate constants (Figs. 7 and 8). The simulations also show that the particular combination of rate constants may result in a transporter that has robust behavior that does not change under different environmental conditions or may switch transport behavior with changing environmental conditions (Fig. 9). This demonstrates the utility of a free-exchange model for systems requiring consideration of “minor” states/transitions and demystifies how substrate and ion can be coupled without restrictions on transporter states and transitions. Nevertheless, free exchange transport is inherently less energetically efficient than stoichiometric transport mechanisms. Here, we consider the potential benefits of a free-exchange model in comparison with a traditional tightly coupled mechanism to begin to understand why some transporters appear to follow each type of mechanism.

Evolution

Evolution selects for coupling tailored to the needs of a particular biological system (Henderson et al., 2019). For example, mammals regulate sugar availability at the organismal level. Consequently, the mammalian glucose transporter (GLUT) family of hexose transporters function as sugar uniporters, with a few tissue-specific exceptions. In contrast, microbes often have to survive in environments with unpredictable nutrient availability. Thus, many microbial homologues of the GLUT transporters appear to be tightly coupled sugar/proton symporters, allowing efficient concentrative sugar import (Walmsley et al., 1998; Wilson-O’Brien et al., 2010; Mueckler and Thorens, 2013). Interestingly, the GLUT homologue GlcP from Staphylococcus epidermidis appears to represent an intermediate case. It employs a free exchange mechanism to perform both proton-coupled sugar symport and sugar uniport depending on external pH (Bazzone et al., 2017). These examples of sugar transporters demonstrate how each have likely evolved to maintain the degree of coupling sufficient for their function in the host organism.

While sugar symporters represent an instance where imperfect coupling is sufficient, loose coupling may be required for multidrug antiporters. Observations of ATP-driven multidrug transporters reveal that futile cycles of ATP-hydrolysis occur between productive drug efflux events and that this may be a trade-off required to achieve multidrug recognition (Al-Shawi et al., 2003; Sharom et al., 1995). Thus, it was predicted that proton-driven multidrug transporters may also sacrifice efficiency to achieve multidrug recognition (Krupka, 1999a,b). The ability to export diverse toxic substrates is critical to ensure survival, and thus would be favored evolutionarily, even at the expense of perfect efficiency. A second potential advantage of a free exchange mechanism involves efficient substrate release. Multidrug transporters recognize and transport chemically diverse molecules, and the binding affinities (and thus off-rates) can range over five orders of magnitude in the case of EmrE (Morrison and Henzler-Wildman, 2014). This presents an issue for tight binding substrates, where slow off-rates may cause tight binding substrates to become “stuck” on the transporter if not for the ability to simultaneously bind proton, which reduces the affinity and enhances the drug off-rate. This idea is supported by Gillespie simulations of EmrE that can resolve the exact pathways for transport (Robinson et al., 2017). This flexibility is not possible in pure-exchange models of transport (Fig. 1) but is possible within the free exchange (Figs. 2 B and 3), and may be a property that is selected for by evolution, by enabling continued efficient transport of toxins out of the cell.

Even for transporters with apparently perfect stoichiometric transport like LacY, single point mutations can produce leak pathways (Forrest et al., 2011; Guan and Kaback, 2006; Varela and Wilson, 1996). The fine-tuning required to achieve stoichiometric transport in naturally occurring transporters may help to explain the relatively loose coupling exhibited by the de novo–designed zinc/proton exchanger Rocker (Joh et al., 2014). However, this suggests that there may be a low barrier for evolution from a tightly to a loosely coupled transport regimen, and thus, tightly coupled systems could become “scaffolding” to design new free-exchange transporters. In fact, this appears to have happened naturally in the small multidrug resistance transporter family (Kermani et al., 2018).

Voltage in the proton-motive force

The proton motive force that drives proton-coupled transport in vivo includes both the proton gradient (ΔpH) and the transmembrane voltage (ΔΨ). Indeed, coupled antiport in EmrE can be driven by either ΔpH or ΔΨ alone or in combination (Rotem and Schuldiner, 2004; Robinson et al., 2018b; Yerushalmi et al., 1995). One limitation of the current free-exchange model is that we consider only ΔpH and not ΔΨ. This choice was based on the availability of experimental measurements of rates for EmrE at different pH values, allowing us to more reliably model the effect of ΔpH. However, since alternating access moves both charged substrates and charged residues of EmrE across the membrane, alternating-access rates are certainly impacted by ΔΨ. In addition, the binding of charged substrates can be affected by voltage as well, if the binding process effectively moves a charged substrate partway down the voltage gradient (Garcia-Celma et al., 2009). Unfortunately, despite the ability of voltage gradients to drive transport, the effect of voltage on individual steps in the transport cycle is unclear (Forrest et al., 2011). New research probing the effect of voltage on all steps in the transport cycle is needed to elucidate how voltage influences the microscopic rate constants required to model transport.

Overall, modeling voltage presents new possibilities to predict transport under more conditions and explain some phenotypes of EmrE. For example, WT EmrE has been seen to move proton and drug in the same direction under voltage alone, suggesting the ability to reverse the flux of traditionally antiported substrates (Robinson et al., 2017). Voltage may also enhance the efficiency of electrogenic coupled transport by further biasing relative flux throughout certain pathways. For example, the negative-inside potential could speed up proton on-rates and slow proton off-rates on the periplasmic side of the membrane, resulting in faster turnover of drug. Overall, more complex phenotypes may occur when both ΔΨ and ΔpH are present. This is the subject of future effort in our laboratory to acquire the necessary kinetic data to expand the model to include voltage in a biologically relevant manner.

Conclusions

Research into free exchange models of transport has broad implication for protein design. By understanding the requirements for different transport regimes, and how to switch between them, researchers can design de novo transporters or manipulate existing transporters to switch behavior under different pH, drug, or other environmental conditions. Overall, our understanding of mechanisms not as a singular pathway but as a superposition of many possible outcomes allows us to further investigate the boundaries of transporters, understand the physiological behavior of these important proteins, and engineer novel systems.

Online supplemental material

The supplemental information includes the matrix form of the differential equations (Eqs. 1 and 2) as well as the Matlab files for the kinetic model. All of the Matlab files for the numerical simulations, including the specific parameters used to produce each of the figures in this paper, are freely available on Github.

Acknowledgments

Merritt C. Maduke served as editor.

We thank Hallie Hanson for assistance writing the Matlab code.

Research reported in this publication was supported by National Institutes of Health grant R01GM0958339 to K A. Henzler-Wildman. G A. Hussey was partially supported by the Hilldale Undergraduate Research Fellowship at University of Wisconsin–Madison and N.E. Thomas by National Institutes of Health grant T32GM007215. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

The authors declare no competing financial interests.

Author contributions: The project was conceived by K A. Henzler-Wildman. Work was carried out by G A. Hussey with advice from K A. Henzler-Wildman and N.E. Thomas. Results were analyzed and manuscript was written by K A. Henzler-Wildman, G A. Hussey, and N.E. Thomas.

References

References
Adam
,
Y.
,
N.
Tayer
,
D.
Rotem
,
G.
Schreiber
, and
S.
Schuldiner
.
2007
.
The fast release of sticky protons: kinetics of substrate binding and proton release in a multidrug transporter
.
Proc. Natl. Acad. Sci. USA.
104
:
17989
17994
.
Al-Shawi
,
M.K.
,
M.K.
Polar
,
H.
Omote
, and
R A.
Figler
.
2003
.
Transition state analysis of the coupling of drug transport to ATP hydrolysis by P-glycoprotein
.
J. Biol. Chem.
278
:
52629
52640
.
Amaral
,
L.
,
A.
Martins
,
G.
Spengler
, and
J.
Molnar
.
2014
.
Efflux pumps of Gram-negative bacteria: what they do, how they do it, with what and how to deal with them
.
Front. Pharmacol.
4
:
168
.
Bakker
,
E.P.
, and
W.E.
Mangerich
.
1981
.
Interconversion of components of the bacterial proton motive force by electrogenic potassium transport
.
J. Bacteriol.
147
:
820
826
.
Bazzone
,
A.
,
A.J.
Zabadne
,
A.
Salisowski
,
M.G.
Madej
, and
K.
Fendler
.
2017
.
A Loose Relationship: Incomplete H+/Sugar Coupling in the MFS Sugar Transporter GlcP
.
Biophys. J.
113
:
2736
2749
.
Bozzi
,
A.T.
,
L B.
Bane
,
C.M.
Zimanyi
, and
R.
Gaudet
.
2019
.
Unique structural features in an Nramp metal transporter impart substrate-specific proton cotransport and a kinetic bias to favor importroton co-transport and voltage dependence enforce unidirectional metal transport in an Nramp transporter
.
J. Gen. Physiol.
doi:
Brill
,
S.
,
O.S.
Falk
, and
S.
Schuldiner
.
2012
.
Transforming a drug/H+ antiporter into a polyamine importer by a single mutation
.
Proc. Natl. Acad. Sci. USA.
109
:
16894
16899
.
Chen
,
Y.-J.
,
O.
Pornillos
,
S.
Lieu
,
C.
Ma
,
A.P.
Chen
, and
G.
Chang
.
2007
.
X-ray structure of EmrE supports dual topology model
.
Proc. Natl. Acad. Sci. USA.
104
:
18999
19004
.
Craciun
,
G.
, and
M.
Feinberg
.
2006
.
Multiple Equilibria in Complex Chemical Reaction Networks: II. The Species-Reaction Graph
.
SIAM J. Appl. Math.
66
:
1321
1338
.
Dohán
,
O.
,
C.
Portulano
,
C.
Basquin
,
A.
Reyna-Neyra
,
L.M.
Amzel
, and
N.
Carrasco
.
2007
.
The Na+/I symporter (NIS) mediates electroneutral active transport of the environmental pollutant perchlorate
.
Proc. Natl. Acad. Sci. USA.
104
:
20250
20255
.
Feng
,
L.
,
E B.
Campbell
, and
R.
MacKinnon
.
2012
.
Molecular mechanism of proton transport in CLC Cl-/H+ exchange transporters
.
Proc. Natl. Acad. Sci. USA.
109
:
11699
11704
.
Fitzgerald
,
G A.
,
C.
Mulligan
, and
J A.
Mindell
.
2017
.
A general method for determining secondary active transporter substrate stoichiometry
.
eLife.
6
:e21016.
Forrest
,
L.R.
,
R.
Krämer
, and
C.
Ziegler
.
2011
.
The structural basis of secondary active transport mechanisms
.
Biochim. Biophys. Acta.
1807
:
167
188
.
Garcia-Celma
,
J.J.
,
I.N.
Smirnova
,
H.R.
Kaback
, and
K.
Fendler
.
2009
.
Electrophysiological characterization of LacY
.
Proc. Natl. Acad. Sci. USA.
106
:
7373
7378
.
Gayen
,
A.
,
M.
Leninger
, and
N.J.
Traaseth
.
2016
.
Protonation of a glutamate residue modulates the dynamics of the drug transporter EmrE
.
Nat. Chem. Biol.
12
:
141
145
.
Guan
,
L.
, and
H.R.
Kaback
.
2006
.
Lessons from lactose permease
.
Annu. Rev. Biophys. Biomol. Struct.
35
:
67
91
.
Henderson
,
R.K.
,
K.
Fendler
, and
B.
Poolman
.
2019
.
Coupling efficiency of secondary active transporters
.
Curr. Opin. Biotechnol.
58
:
62
71
.
Horn
,
F.
, and
R.
Jackson
.
1972
.
General mass action kinetics
.
Arch. Ration. Mech. Anal.
47
:
81
116
.
Iancu
,
C.V.
,
J.
Zamoon
,
S B.
Woo
,
A.
Aleshin
, and
J.Y.
Choe
.
2013
.
Crystal structure of a glucose/H+ symporter and its mechanism of action
.
Proc. Natl. Acad. Sci. USA.
110
:
17862
17867
.
Joh
,
N.H.
,
T.
Wang
,
M.P.
Bhate
,
R.
Acharya
,
Y.
Wu
,
M.
Grabe
,
M.
Hong
,
G.
Grigoryan
, and
W.F.
DeGrado
.
2014
.
De novo design of a transmembrane Zn2+-transporting four-helix bundle
.
Science.
346
:
1520
1524
.
Kermani
,
A A.
,
C B.
Macdonald
,
R.
Gundepudi
, and
R B.
Stockbridge
.
2018
.
Guanidinium export is the primal function of SMR family transporters
.
Proc. Natl. Acad. Sci. USA.
115
:
3060
3065
.
Krupka
,
R.M.
1999
a
.
Limits on the tightness of coupling in active transport
.
J. Membr. Biol.
167
:
35
41
.
Krupka
,
R.M.
1999
b
.
Uncoupled active transport mechanisms accounting for low selectivity in multidrug carriers: P-glycoprotein and SMR antiporters
.
J. Membr. Biol.
172
:
129
143
.
Lolkema
,
J.S.
, and
D.J.
Slotboom
.
2019
.
Models to determine the kinetic mechanisms of ion-coupled transporters
.
J. Gen. Physiol.
151
:
369
380
.
Luo
,
P.
,
S.
Dai
,
J.
Zeng
,
J.
Duan
,
H.
Shi
, and
J.
Wang
.
2018
.
Inward-facing conformation of l-ascorbate transporter suggests an elevator mechanism
.
Cell Discov.
4
:
35
.
Morrison
,
E A.
, and
K A.
Henzler-Wildman
.
2014
.
Transported substrate determines exchange rate in the multidrug resistance transporter EmrE
.
J. Biol. Chem.
289
:
6825
6836
.
Morrison
,
E A.
,
G.T.
DeKoster
,
S.
Dutta
,
R.
Vafabakhsh
,
M.W.
Clarkson
,
A.
Bahl
,
D.
Kern
,
T.
Ha
, and
K A.
Henzler-Wildman
.
2011
.
Antiparallel EmrE exports drugs by exchanging between asymmetric structures
.
Nature.
481
:
45
50
.
Morrison
,
E A.
,
A.E.
Robinson
,
Y.
Liu
, and
K A.
Henzler-Wildman
.
2015
.
Asymmetric protonation of EmrE
.
J. Gen. Physiol.
146
:
445
461
.
Mueckler
,
M.
, and
B.
Thorens
.
2013
.
The SLC2 (GLUT) family of membrane transporters
.
Mol. Aspects Med.
34
:
121
138
.
Nguitragool
,
W.
, and
C.
Miller
.
2006
.
Uncoupling of a CLC Cl-/H+ exchange transporter by polyatomic anions
.
J. Mol. Biol.
362
:
682
690
.
Oh
,
S.
, and
O.
Boudker
.
2018
.
Kinetic mechanism of coupled binding in sodium-aspartate symporter GltPh
.
eLife.
7
:e37291.
Purewal
,
A.S.
1991
.
Nucleotide sequence of the ethidium efflux gene from Escherichia coli
.
FEMS Microbiol. Lett.
66
:
229
231
.
Quistgaard
,
E.M.
,
C.
Löw
,
F.
Guettou
, and
P.
Nordlund
.
2016
.
Understanding transport by the major facilitator superfamily (MFS): structures pave the way
.
Nat. Rev. Mol. Cell Biol.
17
:
123
132
.
Robinson
,
A.E.
,
J.P.
Henderson
, and
K A.
Henzler-Wildman
.
2018
.
A mass spectrometry based transport assay for studying EmrE transport of unlabeled substrates
.
Anal. Biochem.
549
:
130
135
.
Robinson
,
A.E.
,
N.E.
Thomas
,
E A.
Morrison
,
B.M.
Balthazor
, and
K A.
Henzler-Wildman
.
2017
.
New free-exchange model of EmrE transport
.
Proc. Natl. Acad. Sci. USA.
114
:
E10083
E10091
.
Rotem
,
D.
, and
S.
Schuldiner
.
2004
.
EmrE, a multidrug transporter from Escherichia coli, transports monovalent and divalent substrates with the same stoichiometry
.
J. Biol. Chem.
279
:
48787
48793
.
Sharom
,
F.J.
,
X.
Yu
,
J.W.
Chu
, and
C A.
Doige
.
1995
.
Characterization of the ATPase activity of P-glycoprotein from multidrug-resistant Chinese hamster ovary cells
.
Biochem. J.
308
:
381
390
.
Silverman
,
J A.
2002
.
Multidrug-Resistance Transporters
. In
Membrane Transporters as Drug Targets.
G.L.
Amidon
, and
W.
Sadée
, editors.
Springer US
,
Boston, MA
.
353
386
.
Stein
,
W.
1986
.
Transport and Diffusion Across Cell Membranes.
Academic
,
San Diego, CA
.
Thomas
,
N.E.
,
C.
Wu
,
E A.
Morrison
,
A.E.
Robinson
,
J.P.
Werner
, and
K A.
Henzler-Wildman
.
2018
.
The C terminus of the bacterial multidrug transporter EmrE couples drug binding to proton release
.
J. Biol. Chem.
293
:
19137
19147
.
Varela
,
M.F.
, and
T.H.
Wilson
.
1996
.
Molecular biology of the lactose carrier of Escherichia coli
.
Biochim. Biophys. Acta.
1276
:
21
34
.
Walmsley
,
A.R.
,
M.P.
Barrett
,
F.
Bringaud
, and
G.W.
Gould
.
1998
.
Sugar transporters from bacteria, parasites and mammals: structure-activity relationships
.
Trends Biochem. Sci.
23
:
476
481
.
Wilson-O’Brien
,
A.L.
,
N.
Patron
, and
S.
Rogers
.
2010
.
Evolutionary ancestry and novel functions of the mammalian glucose transporter (GLUT) family
.
BMC Evol. Biol.
10
:
152
.
Yazaki
,
K.
,
A.
Sugiyama
,
M.
Morita
, and
N.
Shitan
.
2008
.
Secondary transport as an efficient membrane transport mechanism for plant secondary metabolites
.
Phytochem. Rev.
7
:
513
524
.
Yerushalmi
,
H.
,
M.
Lebendiker
, and
S.
Schuldiner
.
1995
.
EmrE, an Escherichia coli 12-kDa multidrug transporter, exchanges toxic cations and H+ and is soluble in organic solvents
.
J. Biol. Chem.
270
:
6856
6863
.
Zhang
,
Y.
,
Y.
Zhang
,
K.
Sun
,
Z.
Meng
, and
L.
Chen
.
2019
.
The SLC transporter in nutrient and metabolic sensing, regulation, and drug development
.
J. Mol. Cell Biol.
11
:
1
13
.
This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms/). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 4.0 International license, as described at https://creativecommons.org/licenses/by-nc-sa/4.0/).

Supplementary data