Modeling Bacterial Flagellar Motor With New Structure Information: Rotational Dynamics of Two Interacting Protein Nano-Rings

In this article, we develop a mathematical model for the rotary bacterial flagellar motor (BFM) based on the recently discovered structure of the stator complex (MotA5MotB2). The structure suggested that the stator also rotates. The BFM is modeled as two rotating nano-rings that interact with each other. Specifically, translocation of protons through the stator complex drives rotation of the MotA pentamer ring, which in turn drives rotation of the FliG ring in the rotor via interactions between the MotA ring of the stator and the FliG ring of the rotor. Preliminary results from the structure-informed model are consistent with the observed torque-speed relation. More importantly, the model predicts distinctive rotor and stator dynamics and their load dependence, which may be tested by future experiments. Possible approaches to verify and improve the model to further understand the molecular mechanism for torque generation in BFM are also discussed.


INTRODUCTION
The swimming motion of many bacteria is powered by the rotary bacterial flagellar motor (BFM). As a canonical example, the motors in E. coli are fueled by electrochemical potential difference of protons (H + ) across the cytoplasmic membrane, i.e., the proton motive force (PMF) (Larsen et al., 1974;Hirota et al., 1981;Berg, 2003). The BFM can rotate bidirectionally. Without binding to the phosphorylated response regulator CheY protein (CheY-P), the BFM rotates in a counterclockwise (CCW, viewed from the filament to the motor) direction. There are ∼ 4 − 6 BFMs in an E. coli cell. The CCW rotation of the BFMs causes their flagellar filaments to form a coherent bundle and propels the swimming cells in a roughly straight line that is called a "run." However, the binding of CheY-P to BFM causes it to rotate in the clockwise (CW) direction, which disrupts the filament bundle. As a result, the rod-shaped E. coli cell changes its orientation randomly without translational movement in an event called a "tumble." The bacterial flagellar motor consists of a rotor and a number of stators that can vary from 0 to 11 depending on gene expression and mechanical load (Blair and Berg, 1988;Lele et al., 2013;Tipping et al., 2013;Nord et al., 2017;Wadhwa et al., 2019Wadhwa et al., , 2021. The rotor consists of FliF, FliG, FliM, and FliN. FliF forms the MS-ring embedded within the cytoplasmic membrane (Asai et al., 1997;Sato and Homma, 2000;Yorimitsu et al., 2004;Morimoto and Minamino, 2014;Baker et al., 2016). The cytoplasmic face of the MS-ring is attached to the C-ring, which comprises the FliG, FliM, and FliN proteins. The diameter of the C-ring is ∼ 45nm. The C-ring is responsible for coupling the rotor to the stator units and for switching rotational directions. The membrane-bound stators, which consist of the MotA and MotB proteins, are powered by the PMF to drive the rotation of the rotor via direct interaction between MotA and FliG (Chun and Parkinson, 1988;Blair and Berg, 1990;Kojima and Blair, 2004;Roujeinikova, 2008).
The molecular mechanism of torque generation by the stator units has been suggested by two recent Cryo-EM studies of the stator structure (Deme et al., 2020;Santiveri et al., 2020). Both studies report that a stator unit consists of a MotB dimer surrounded by a MotA pentamer ring with a diameter of ∼ 5 − 7.5nm. The position of the MotB dimer, which is bound to the peptidoglycan layer in the rigid cell wall, remains fixed, while the MotA pentamer forms a ring that can rotate around the MotB dimer. Each MotB monomer progressively engages with each MotA monomer to form an ion channel through which one proton can pass to drive the rotation of the MotA pentamer ring.
Brownian ratchet models have attracted physicists' attention for modeling molecular motors since Richard Feynman popularized the idea more than half a century ago (Feynman et al., 1966;Parrondo et al., 1998;Astumian, 2010;Golubeva et al., 2012). Among all variants of the ratchet models, only isothermal chemical ratchets are relevant for biological motors (Jülicher et al., 1997;Parrondo and de Cisneros, 2002). A minimal model for the rotational ratchet model can be found in Xing et al. (2006), Tu (2009), andMeacci et al. (2011). This model has two types of coordinates: the "physical" coordinate that describes the motion of the motor in physical space (angles for rotary motor); and the "chemical" coordinate that describes the different conformational states of the motor. For the BFM, the motor moves in a periodic potential in which torque, which is the gradient of the potential, is generated as the motor slide down from a position of high potential energy toward its equilibrium position, which has the lowest potential energy. The passage of an ion triggers a conformational change (stepping) of the motor. Because the landscape of potential energy is different for different conformational states of the motor, the conformational change driven by the PMF changes the potential energy landscape and can bring the motor to a position of high potential energy in the new potential energy landscape. The newly gained potential energy continues to drive the physical rotation of the motor. This continuous process drives the system toward a sequence of new equilibrium positions and gives rise to a directed stepwise rotation (Sowa et al., 2005).
One of the most important characteristics of a rotary motor is its torque-speed dependence (or force-speed dependence for linear motors). The measured torque-speed curve for the CCW BFM has a downward concave shape with a roughly constant high torque at low to medium speeds and a rapid and nearly linear decrease of torque at high speeds (Chen and Berg, 2000;Lo et al., 2013). This concave torque-speed curve has the advantage of generating high power output (or, equivalently, a high torque for a given speed) in a wide range of physiologically relevant loads. The minimal ratchet model can reproduce the observed torque-speed curve when certain general constraints for the rotor-stator interaction potential and the proton-assisted transition rate between different conformational states are applied (Tu and Cao, 2018).
Before the recent availability of an accurate picture of the structure of the stator (Deme et al., 2020;Santiveri et al., 2020), previous structure-based modeling studies assumed that the stator unit has a MotA 4 MotB 2 composition (Braun et al., 2004;Kim et al., 2008). For example, it was proposed (Mandadapu et al., 2015;Nirody et al., 2016) that the power stroke for each rotational direction of the rotor (CCW or CW) is delivered by two of the four MotA subunits through a steric force generated by a "kink and swivel" conformational change of MotA induced by the binding of a proton to a corresponding MotB. The newly discovered structure of the stator complex has the different stoichiometry of MotA 5 MotB 2 (Deme et al., 2020;Santiveri et al., 2020). The structure of the stator complex immediately led to an important new insight into the molecular mechanism by which the BFM operates: the MotA ring in the stator unit rotates. Consequently, all of the five MotA subunits are involved in torque generation independent of the rotational direction of the rotor, and the power stroke is generated directly by the rotational motion of the MotA subunits themselves instead of their internal conformational changes.
In terms of modeling, the new insight gained from the stator complex structure means that the stator unit must be described by its physical state in terms of the relative angle between the MotA pentamer ring and the MotB dimer in addition to its chemical (conformational) states. In the following section, we present an updated model for the BFM based on the general modeling framework described before (Tu and Cao, 2018), but now in the context of the new structure of the stator complex. Some preliminary results predicted by the model will also be described.

A STRUCTURE-INFORMED MODEL AND SOME PRELIMINARY RESULTS
Our new model describes the dynamics of the two coupled rotating nano-rings (the MotA ring and the FliG ring) as shown in Figures 1A,B. The stator unit is the active part of the motor. It can generate torque as a rotational ratchet powered by the PMF. Rotation of the MotA ring drives the rotation of the passive FliG ring through a coupling potential V R (θ S , θ R ), where θ S and θ R are the rotational angles of the stator and the rotor, which have a periodicity of θ S = 2π/5, θ R = 2π/26, respectively. It is assumed that there are 26 FliG proteins in the C-ring (Thomas et al., 2006). The rotational motion of the rotor is described by a Langevin equation for the rotor angle θ R : where ξ R is the load on the rotor, the first term on the right hand side of Equation 1 is the torque exerted on the rotor from the stator via the MotA-FliG interaction, and η R (t) is the rotor thermal noise: shows the two interaction potentials between the MotA pentamer ring and the MotB dimer: V 1 (θ S ) from MotB 1 (blue) and V 2 (θ S ) from MotB 2 (red), which have the same shape but are shifted from each other by half a period θ S /2 = 2π/10. The blue and red energy landscapes are shifted vertically only for clarity. P and C stand for the two conformational states of MotB in which the proton acceptor residue faces the periplasm and the cytoplasm, respectively, and I stands for an intermediate state between P-state and C-state. The plus (or minus) superscript indicates the proton-bound (or proton-unbound) state. Transition (stepping) from one MotA-MotB interaction potential to another interaction potential occurs with a probability rate k m near the potential minimum (shaded region in the left panel), which is driven by the proton discharge-recharge chain reaction P + → P − → C − → C + . For brevity, only the initial and final states (P + and C + ) are shown, and the intermediate discharged states C − and P − are omitted (see text for the full description). The rotor energy landscape (right panel) shows the interaction potential (V R ) between MotA and FliG, which depends only on the scaled angle difference ϕ = θ R − kθ S . For simplicity, we used a piece-wise linear form for V 1,2 and V R . To the left of the P → I plateau in V 1,2 , a barrier is added to prevent back flow. The four stages of stator dynamics during one full power stroke as described in Equation 2  the Boltzmann constant and T the temperature. The coupling potential V R (θ S , θ R ) should be periodic in both θ S and θ R . To satisfy this dual periodicity, we write , which is the gear ratio between the two rotating rings in E. coli.
We model the stator unit as a rotational ratchet driven by the PMF. Based on the new structural information, the MotA pentamer ring rotates around the fixed MotB dimer. Given the incommensurate stoichiometry (5 : 2) within the stator unit, the two MotB monomers engage with MotA monomers alternately -when one MotB is engaged directly with one MotA monomer the other MotB is positioned between two MotA monomers. Consequently, the two MotB monomers generate two different energy potentials shifted by half a period ( θ S /2). To account for this effect, the overall stator potential can be written as the sum of two identical periodic potentials shifted by an angle difference of Figure 1B.
During one power stroke, the MotA pentamer ring rotates one period (2π/5). Each full power stroke consists of two half strokes, each of which is powered by a proton passing through one of the two ion channels alternately. During each half stroke, the proton acceptor residue (D32 for E. coli) on the corresponding MotB monomer becomes protonated, which drives rotation of the MotA ring. At the same time, the protonated MotB undergoes a conformational change in which the proton acceptor residue changes from a periplasm-facing position to a cytoplasm-facing position, from which the proton can be released to the cytoplasm. To describe this process, we define three conformational states for each MotB: a periplasm-facing state called the P-state, a cytoplasm-facing state called the C-state, and an intermediate state between the P-state and C-state, which we call the I-state.
Each MotB can be either protonated or unprotonated (indicated by a plus or minus sign). A proton can bind with MotB when it is in the P-state; and a bound proton can enter the cytoplasm when the MotB is in the C-state.
As suggested by the recent Cryo-EM study (Santiveri et al., 2020), different conformational states occur at different relative positions between MotB and the closest MotA monomer, or, equivalently, at different locations within the MotA-MotB landscape of interaction energy potential, as shown in Figure 1B (left). In particular, the P-state occurs near the peak potential and the C-state occurs near the minimum. For simplicity, we used a simple piece-wise linear potential (V 1,2 ) that is flat over half of the period and linear in the other half, as shown in Figure 1B (left). Using this simple representation of the internal chemical states of the stator, evolution of the internal states of the two MotB monomers in a stator unit during one full rotational period (2π/5) can be described by the reaction cycle: Proton discharge and recharge on MotB 1 Proton discharge and recharge on MotB 2 I + 1 , (2) where the subscripts (1,2) indicate the two MotB monomers. For each MotB, torque is generated during the transitions P + → I + → C + , which drives the physical rotation of the MotA ring.
After reaching the C + -state near the potential energy minimum, MotB is rejuvenated by the release of the proton to the cytoplasm followed by a conformational change from the C-state to the Pstate and subsequent addition of a new proton to the P-state: C + → C − → P − → P + . This chain reaction returns MotB to its high-energy protonated P-state (P + ) to continue driving the rotation of the stator. It is clear from the reaction cycle (Equation 2) that one full power stroke is powered by translocation of two protons (or 2PMF), each "consumed" by one MotB. The four stages of the alternating physical rotations and proton discharge-charge reactions during one full power stroke as described in Equation 2 are illustrated in Figures 1C-F, respectively. In our model, the rotation dynamics of the MotA ring can also be described by a Langevin equation for the stator angle θ S : where ξ S is the load (viscosity) of the stator, the first and the second terms on the right hand side of Equation 4 are the torque exerted on the stator from the rotor via the MotA-FliG interaction and the torque generated by the MotA-MotB interaction, respectively, and η S (t) is the stator thermal noise with < η S (t)η S (t ′ ) >= 2k B Tξ S δ(t − t ′ ). The shape of V S (θ S ) is given in Figure 1B. To prevent back flow at the maximum potential energy, we introduce an additional energy barrier which extends to the left of the P → I plateau.
The key chemical reaction that powers rotation of the stator is the proton discharge-recharge chain reaction (C + → C − → P − → P + ), which leads to the transition (stepping) from one potential landscape, in which the system has a low (near minimum) potential energy, to another potential landscape where the system has a high (near maximum) potential energy. For simplicity, this chain reaction is modeled here by an overall transition rate k m (θ S ) that depends on the MotA-MotB relative angle θ S (see Figure 1B for an illustration of k m , the blue arrow in the left panel). Specifically, Prob (transition from low to high potential energy during In summary, our model considers the coupled rotations of both the active stator and the passive rotor as well as the key chemical reactions that power the rotation of the stator. The dynamics of the BFM are fully described by Equations 1, 3 and 4, which depend on three biologically meaningful functions: two interaction energy (potential) functions -V R (θ S − kθ R ) between the FliG-ring of the rotor and the MotA ring, V S (θ S ) between the MotB dimer and the MotA ring, and the chemical transition rate function k m (θ S ). Next, we present some preliminary results based on simulations of the model with simple choices for these functions to demonstrate the general utility of the model.
In our simulations, we set the stator load ξ S to be a small constant and change the rotor load ξ R to obtain the torque and speed for different values of ξ R . The motor dynamics depend on the relative strength between the MotA-FliG coupling interaction potential V R and the torque-generating potential V S . In Figure 2, we present the results for a case in which these two potentials are comparable. As shown in Figure 2A, the torque-speed curve shows a downward concave shape consistent with experiments. In addition to the torque-speed relationship, our model allows us to investigate the rotational dynamics of both the rotor and the stator. In Figure 2B, we plot the time trajectories of the stator (dashed lines) and rotor (solid lines) angles for different loads. At a low load, the two trajectories are close to each other, whereas, at a high load, the trajectory of the rotor is almost flat, even though the stator still rotates at a finite speed. As shown in Figure 2C, as load increases, the rotor's speed quickly decreases to zero, while the stator's speed approaches a constant speed (reduced but still finite). Our results indicate that the rotations of the stator and the rotor are in-sync at low load: when the stator rotates at an angle of θ , the rotor rotates at approximately k θ , i.e., there is no slippage between the two rotating rings. However, slippage between the two rings increases with the load. Near stall, the rotor speed vanishes while the stator still rotates at a finite speed. The existence of the stator-rotor slippage depends on the strength of the stator-rotor (MotA-FliG) coupling (V R ) relative to that of the MotA-MotB coupling (V S ). When we used a V R that is much stronger than V S , slippage disappears. Here, we choose to show a case in which V R and V S are comparable to highlight the possibility of slippage between the two rotating protein nanorings.
To investigate the relative motion (rotation) of the stator and rotor in details, we plot the trajectories of θ R (t) vs. θ S (t) for different loads in Figure 2D, where the dashed black line indicates the tight coupling (no slippage) limit when θ R = kθ S (dotted black line). At low load (the green line), the trajectory is almost a straight line parallel to the tight coupling limit, while at high load (the red line), the trajectory wiggles closer to the horizontal direction. We can categorize the relative statorrotor motion by coarse-graining the trajectories in the θ S − θ R space. We found three types of relative motion, which are shown in the zoom-in plots in Figure 2E. In type I movement, the rotor moves with the stator without slippage; in type II movement, the rotor-stator system follows a biased random walk in the (θ S , θ R ) space with net motion lower for both the stator and the rotor; in type III movement, the stator rotates at a finite speed but the rotor has almost no net motion, which indicates a high degree of slippage between the rotor and the stator. Type III movement can be considered as the extreme case of type II movement. From the results of our simulation, we find that as load increases, the motor spends more time in the type II&III movements, i.e., slipping. This observation is further confirmed by the distributions of the scaled angle difference between the rotor and stator, i.e., ϕ = θ R − kθ S , as shown in Figure 2F. At low load, the distribution is concentrated around ϕ = 0, indicating the absence of slippage. As load increases, the distribution of ϕ becomes more populated in the region with ϕ < 0, which shows that the stator leaves the rotor behind more often, i.e., the motor spends more time slipping.
The most unique feature of our model is its ability to investigate the rotational motion of the stator in addition to the rotational motion of the rotor, which was the main focus of previous modeling studies. The model predictions, such as the dependence of the stator speed on the rotor load ( Figure 2C) and the distributions of the relative rotor-stator angle (Figure 2F), can be tested by future experiments that directly measure the rotational motion of the stator in addition to the existing measurements of the rotor motion. These experiments together with our model can also be used to infer the interaction potentials (V S and V R ) as well as the proton-assisted transition rate (k m ).

DISCUSSION
In this article, we propose a new model for BFM based on a general modeling framework and the recently discovered molecular structure of the stator complex. The most important update in our model is that the stator is now modeled as an active rotating nano-machine that consumes PMF and drives the rotation of the rotor. The model can generate the characteristic torque-speed curve observed in experiments. More importantly, our model can be used to investigate and predict the dynamics of both the rotor and the stator, which can be tested in experiments.
It is generally believed that proton flux and motor rotation are "tightly coupled" in the BFM, which implies that there is no proton flux at stall and the efficiency of the BFM is close to 1 near stall (Berg, 2003). Here, we find that if the MotA-FliG interaction is comparable to (or weaker than) the MotA-MotB interaction, there is slippage between the rotor and the stator. The probability of slippage increases as the load increases, which leads to a motor efficiency that is less than 1 near stall. On the other hand, we also find that slippage can be suppressed by increasing the rotor-stator interaction strength. Thus, our study poses a simple but important question. Does the BFM slip? Ultimately, this question needs to be answered experimentally by measuring the rotation of the stator (MotA-ring) at different loads, especially at the high load where the probability of slippage is expected to be the highest if it exists. Alternatively, one may be able to look for signatures of slipping by measuring fluctuations of the rotor's instantaneous speed at a high load. In the absence of slipping, the rotor should run smoothly at a high load with the dominant fluctuations coming from thermal fluctuations (Meacci and Tu, 2009). However, if there is slipping, then fluctuations in the rotor speed can be much larger with additional contributions coming from the slipping events.
Our model can be extended to cases with multiple stators. At low load, the stators and rotor rotate in sync without slippage, so the rotor's maximum speed is determined by the maximum speed of the individual stator units independent of the number of stators. However, at high load (e.g., near stall), the effective duty ratio of the stator can be far below 1 due to the higher probability of slippage and/or the relatively long time spent at the energy potential minima waiting for the conformational change (Meacci and Tu, 2009). Note that the effective duty ratio during slipping is zero although the stator is not disengaged with the rotor. Assuming no direct coordination among stators, there should be little overlap among the duty cycles of each stator. Thus, the maximum torque at stall is simply the sum of torques generated by all the stators. The dependence of maximum speed and maximum torque on the number of stators has been reported before (Blair and Berg, 1988;Yuan and Berg, 2008), and our model may provide a natural explanation for these experimental observations.
The rotation of the stator units is also essential for understanding switching between CW and CCW rotation of the rotor. As pointed out in Santiveri et al. (2020) and illustrated in Figure 1A here, the MotA ring always rotates clockwise around the MotB dimer, and the switch between the CW and CCW rotation of the rotor ring is determined by conformational changes of the C-ring. In the CCW mode, the C-ring adopts a compact conformation, and the stator units contact FliG subunits on the outside of the C-ring; thus, the C-ring rotates CCW. Upon binding of CheY-P to FliM in the C-ring, the C-ring expands, and the stator units are now in contact with FliG subunits on the inside of the C-ring, so the C-ring rotates clockwise. In this manner, the unidirectional rotation of the stator drives the bidirectional rotation of the motor. In our modeling framework, the coupling between rotor and stator is modeled by V R (θ R , θ S ). We hypothesize that the CW and CCW rotation modes correspond to different forms of the interaction potential V R , which can lead to the different torque-speed curves for the BFM in CCW and CW rotation observed experimentally (Yuan et al., 2010).
The BFM, like all other biological motors, converts chemical energy (e.g. from ATP hydrolysis or ion translocations) into mechanical work. A fundamental question is whether there are thermodynamic bounds to the power generation and energy efficiency for these highly non-equilibrium molecular engines (Parmeggiani et al., 1999;Parrondo and de Cisneros, 2002;Astumian, 2010). A related question is what microscopic design principles allow a molecular motor to approach these bounds under realistic constraints. In our model, the motor design space is spanned by three functions: V S (θ S ), V R (ϕ), and k m (θ S ), which determine the motor dynamics. For example, as we showed in this paper, depending on the relative strength of V R and V S , there could be slippage between the stator and the rotor, which strongly affects the efficiency of the motor. Each of these functions contains information about the molecular details of the motor, which, in principle, can be estimated from or constrained by the structural data. A systematical exploration of the motor design space will lead to a better understanding of the fundamental limits and possible design principles for optimal motor performance.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
YT originated the study. YC and TL performed all simulations. YC and YT prepared the manuscript. All authors read and approved the final manuscript.

FUNDING
The work by YT is supported by a NIH grant (R35GM131734).

ACKNOWLEDGMENTS
YT would like to express deepest thanks to Howard Berg for introducing the bacterial flagellar motor-a biotechnological marvel (Berg, 2003) to him more than a decade ago and for Howard's constant enlightenment and encouragement on the topic till his untimely passing at the end of 2021. We dedicate this work in the memory of Howard-it was Howard and members of the Berg lab who first told us about the new stator structure (Deme et al., 2020;Santiveri et al., 2020), who launched this study. YT would also like to thank Nicholas Taylor for an insightful discussion on details of the stator structure and proton translocation process.

SUPPLEMENTARY MATERIAL
Piece-wise linear forms for V R and V 1,2 are used. V 1,2 : k B T = 4.11pN × nm, V 0S = 10k B T, τ + S = 2V 0S / θ S , τ − S = 4V 0S / θ S , where τ ± S is the positive and negative torque on stator potential V 1,2 , respectively, and the extended potential to the left of P → I has the same slope as τ + S . V R : a symmetric form is used with V 0R = 10k B T. Load on the stator ξ S = 0.01(pN nm rad −1 s). The transition rate is k m = 10 5 s −1 for θ S ≥ 2π/5 and 0 otherwise. Equations 1 and 3 are solved with discrete time step dt = 5 × 10 −6 , and in the format of θ i+1 R,S = θ i R,S + dt(· · · )/ξ R,S + 2k B Tdt/ξ R,S N(0, 1), where (· · · ) is the drift term, and N(0, 1) is a random number with standard normal distribution. During each step, the transition from low to high energy potential happens when k m (θ )dt > U(0, 1), where U(0, 1) is a random number with uniform distribution on (0,1).