Fast crystallization of rotating membrane proteins

TL;DR

This paper demonstrates that actively rotating membrane proteins can self-assemble into stable, rotating hexagonal lattices driven by hydrodynamic and repulsive interactions, revealing a novel mechanism for membrane protein crystallization.

Contribution

The study introduces a new understanding of how rotational hydrodynamic interactions facilitate the crystallization of membrane proteins into stable lattices.

Findings

Rotating membrane proteins can form hexagonal lattices.

The crystal rotation speed depends on protein density and lattice size.

Rotational interactions act as a bounded, stabilizing 'temperature' for order.

Abstract

We examine the interactions between actively rotating proteins moving in a membrane. Experimental evidence suggests that such rotor proteins, like the ATP synthases of the inner mitochondrial membrane, can arrange themselves into lattices. We show that crystallization is possible through a combination of hydrodynamic and repulsive interactions between the rotor proteins. In particular, hydrodynamic interactions induce rotational motion of the rotor protein assembly that, in the presence of repulsion, drives the system into a hexagonal lattice. The entire crystal rotates with an angular velocity which increases with motor density and decreases with lattice diameter - larger and sparser arrays rotate at a slower pace. The rotational interactions allow ensembles of proteins to sample configurations and reach an ordered steady state, which are inaccessible to the quenched nonrotational system. Rotational interactions thus act as a sort of temperature that removes disorder, except that actual thermal diffusion leads to expansion and loss of order. In contrast, the rotational interactions are bounded in space. Hence, once an ordered state is reached, it is maintained at all times.

Equations

$$0={\eta }_{2D}{\nabla }^2\mathbf{v}+{\eta }_{3D}{\left[\frac{\partial \mathbf{u}}{\partial z}\right]}_{z=0}+\tau {\nabla }^{\perp }\delta (\mathbf{r}),$$

$$\tilde{v}(q)=\frac{\tau }{{\eta }_{2D}}{\nabla }^{\perp }\tilde{\Psi };\tilde{\Psi }=\frac{1}{q(q+{\lambda }^{-1})}.$$

$$\mathbf{v}(\mathbf{r})$$

$$\Psi (\mathbf{r})$$

$$\frac{d{z}_k}{dt}=i\sum _j\frac{z_k-z_j}{|z_k-z_j{|}^3},$$

$$\Omega \stackrel{N\to \infty }{\to }\frac{4}{D^3}\sum _{j=0}^{n_r}\frac{\cos \left(\frac{\pi n}{3n_r}\right)}{n_r^2}=\frac{6\sqrt{3}}{\pi D^3{n}_r}$$

Full text

A biological membrane is not a homogeneous fluid; it includes around thirty percent transmembrane proteins, making the membrane a quasi-2D suspension Oppenheimer and Diamant (2009). These proteins play various roles and exhibit myriad behaviors. Some proteins passively diffuse in the plane of the membranes Levine et al. (2004), while others are actively driven Gov (2004); Hosaka et al. (2017). In this Letter, we examine the flows and patterns created by one very important transmembrane protein — ATP synthase, a protein that creates the cellular fuel ATP (adenosine triphosphate) by phosphorylation of ADP (adenosine diphosphate), and which rotates as part of its function Davies et al. (2012). ATP synthase has two counter-rotating structures, the ${\rm F_{O}}$ subunit which sits primarily in the membrane, and the ${\rm F_{1}}$ subunit which protrudes into the outer fluid; see Fig. 1a. The rotation rates are proportional to the flux of protons through the protein (and thence across the membrane) and range between $10^{2-3}$ Hz Ueno et al. (2005); Martin et al. (2014). Such rotor proteins are a self-driven active-matter system Nguyen et al. (2014); Yeo et al. (2015); Soni et al. (2018). Since no external torque acts upon the protein its induced flow resembles a torque-dipole to a leading order Lenz et al. (2004).

Electron cryomicroscopy of ATP synthase in lipid vesicles shows that rotor proteins can form lattices in the flat regions of the vesicle Jiko et al. (2015). By studying infinite systems of model rotor proteins, it was suggested earlier by Lenz et al. Lenz et al. (2003, 2004) that a hexagonal lattice is a neutrally stable state. For finite systems we show that adding repulsive interactions leads to a lattice state even when starting from random initial conditions. Figure 1b–d shows the central result of this work — snapshots of three different simulations of rotor protein assemblies, which interact via different combinations of repulsion and rotation: panel b (marked \scalerelS), only hydrodynamic interactions; panel c (marked \scalerelS), only repulsive interactions; and panel d (marked \scalerelS), both hydrodynamic and repulsive interactions. The insets show the structure factor of each system, given by $S({\bf q})=\sum_{i,j}e^{-i{\bf q}\cdot({\bf r}_{i}-{\bf r}_{j})}$ Chaikin and Lubensky (1995). The six distinct peaks in the structure factor of System \scalerelS indicate the presence of global hexagonal order, versus the lack of distinct global features in the structure factors of Systems \scalerelS and \scalerelS.

The hydrodynamics model and its properties: To understand the dynamics of a rotor protein assembly, let us describe a simpler starting point — a membrane rotor, modeled as a single disk of radius $a$ rotating due to an external torque $\tau$, and sitting in a flat membrane whose velocity is ${\bf v}$. We assume membrane incompressibility ($\nabla\cdot{\bf v}=0$), and negligible inertia (small Reynolds number). Under these assumptions, momentum conservation in the membrane reads,

[TABLE]

where ${\bf u}$ is the 3D flow in the outer fluids, and ${\nabla^{\perp}}=(-\partial/\partial y,\partial/\partial x)$. The second term on the right hand side is the jump in shear stress from the outer fluids, and the third term is the force due to a point torque. There is no pressure contribution for purely rotational motion, or a superposition of such flows. The outer fluids obey the 3D Stokes equations with the boundary condition ${\bf u}^{\pm}|_{z=0}={\bf v}$. It is easy to solve this coupled system for the 2D stream-function $\Psi$ of the velocity ${\bf v}$ using a 2D Fourier transform ($\tilde{F}({\bf q})=\int\int F({\bf r})e^{i{\bf q}\cdot{\bf r}}d^{2}r$), giving:

[TABLE]

In real space $\Psi({\bf r})=1/4(H_{0}({r/\lambda})-Y_{0}({r/\lambda}))$, where $H_{0}$ and $Y_{0}$ denote the order zero Struve function and Bessel function of the second kind, respectively.

To derive the stream function for a rotor protein, a second counter-rotating disk is placed a distance $l$ away in the outer fluid, such that the total torque on the protein is zero (see Figure 1A). A similar, albeit lengthier derivation, assuming $r\gg l$ (see Supplementary Information and Ref. Lenz et al. (2003)), gives an almost identical result for the stream function but with one additional term,

[TABLE]

where $\Gamma=\frac{2\tau l}{{\eta_{2D}}\lambda}$. We can now consider the scaling of Eq. (3) in the various limits. For small distances compared to the Saffman-Delbrück (SD) length, $r\ll\lambda$, the stream function satisfies $\Psi\sim 1/r$. Therefore $v\sim{\nabla^{\perp}}\Psi\sim 1/r^{2}$. Qualitatively, this can be understood as follows: in a 2D fluid, a Stokeslet (the flow due to a point force $\delta({\bf r})$) scales as $\log r$; a rotlet (the flow due to a point torque ${\nabla^{\perp}}\delta({\bf r})$) thus scales as $v\sim 1/r$. The flow due to a torque dipole scales as yet another derivative, leading to $v\sim 1/r^{2}$. In the opposite limit, $r\gg\lambda$, the 3D fluid dominates. The flow due to a point force at large distances scales as $1/r$, we would expect the flow due to a rotating protein to scale as $1/r^{3}$, but this term cancels due to symmetry at the $z=0$ plane, and so we have $v\sim 1/r^{4}$.

For more than one rotor protein, the stream function can be generalised to the Hamiltonian $\mathcal{H}=\sum_{i\neq j}\Gamma_{i}\Gamma_{j}\Psi(|r_{i}-r_{j}|)$, where $\Gamma_{i}=\frac{2\tau_{i}l}{{\eta_{2D}}\lambda}$ is the strength of the $i$th torque dipole. The velocity is given by ${\bf v}_{i}=(1/\Gamma_{i})\nabla_{i}^{\perp}\mathcal{H}$ Lushi and Vlahovska (2015). Note that $\mathcal{H}$ is a Hamiltonian in the 2D coordinates, $x$ and $y$, (unlike in classical mechanics where it is a function of position and conjugate momentum), and so phase-space corresponds to the positions of the rotors. From Noether’s theorem Noether (1918), symmetries of the Hamiltonian correspond to conservation laws. In our case there is conservation of the Hamiltonian itself, and of the first- and second-moments (from time, translational and rotational invariance respectively).

We take all the “circulations”, $\Gamma_{i}$, to be equal, $\Gamma_{i}=\Gamma$, and hence same-signed111For ATP synthase, angular velocity of the proteins is a function of proton concentration. Since protons are seven orders of magnitude smaller than the proteins, we can assume that proton concentration equilibrates on a much faster time-scale.. Conservation of the second-moment then simplifies to $\mu_{2}=\sum_{i,j}|r_{i}-r_{j}|^{2}={\rm Const}$. As a result of this invariance the distance between rotor proteins cannot diverge to infinity, and by invariance of $\mathcal{H}$ the distance cannot collapse to zero. In general, rotor protein ensembles stay bounded in an area not much different than their initial area, and maintain a finite distance between each other. A crude upper bound for the radius of the rotating system can be derived as follows. For a random initial condition, the second moment is proportional to the initial area and the number of particles squared, i.e. $\mu\sim r_{0}^{2}N^{2}$. The maximal radius is bounded by a configuration in which all particles but one are at the origin. The remaining particle will then have the maximal possible distance from the origin, given by $r_{\rm max}\leq\sqrt{\mu_{2}/N}\propto r_{0}\sqrt{N}$. This bound can be improved by incorporating conservation of the Hamiltonian, see Supplementary Information.

Surprisingly, the dynamics of such a system is very similar to that of point singularities in 2D such as the ideal vortices of a 2D Euler fluid, or those of the quasi-geostrophic (QG) equations which arise in modeling atmospheric flows Falkovich (2009); Conti and Badin (2017). In all such systems two singularities will orbit around each other, and for four or more the system becomes non-integrable and the dynamics can be chaotic Aref and Pomphrey (1982) (see SI for a video of typical rotor protein dynamics). Indeed, the near-field interactions for rotor proteins are the same as for QG vortices.

In what follows we consider two aspects of membrane proteins that are not captured by pure hydrodynamic interactions. First, membrane proteins are physical objects of finite size with possible interactions with other proteins. Second, due to their small size ($\sim 10$ nm) thermal noise might play a significant role in their dynamics. We show that adding any type of repulsion between the rotor proteins can result in the formation of crystals; sufficiently high temperature can destroy that order.

Repulsive interactions drive crystalization. To start, we follow Lenz et al. (2003) who noted that for a system with only hydrodynamic interactions (System \scalerelS), an infinite hexagonal array of rotor proteins is a steady-state of the system, for any lattice scale $d$. Given the system’s Hamiltonian structure, this array is, at most, neutrally stable. This means that if the system does not start in an ordered array it will not reach that state. For more details see the SI. Adding repulsive interactions (System \scalerelS) can turn this neutrally stable fixed point into a stable one. For example, a linear perturbation, $\delta r(t)$, on a single rotor protein must have the form, $d\delta{\bf r}/dt=-\alpha\delta{\bf r}_{\perp}-\beta\delta{\bf r}$, where $\alpha$ ($\beta$) is coming from the rotational (repulsive) interaction. The eigenvalues of this system are $\lambda=\pm i\alpha-\beta$. The additional negative component is a necessary condition for linear stability of the hexagonal configuration. From simple symmetry arguments it is clear that an infinite or a confined system with long ranged repulsive interactions must form a hexagonal lattice, also known as a Wigner crystal Wigner (1934).

However, this is not the case for a finite system with only short-ranged repulsive interactions, as can be seen from Fig. 1c. Surprisingly, experiments by Jiko et al. Jiko et al. (2015) still show an ordered formation even for finite ensembles of proteins. We consequently examined the dynamics of a finite, unconfined system of rotor proteins with short-ranged repulsion with an interaction distance $r_{s}$, using either exponential (e.g. electrostatically screened interactions) $U=U_{0}e^{-r/r_{s}}$, or soft, harmonic-like repulsion of the form $U=U_{0}(r-r_{s})^{2}$ up to an inter-particle distance of $r_{s}$, and zero otherwise.

Starting from random initial conditions and letting the system evolve, we discovered that rotation promotes rapid organization into crystals even in cases where the repulsion alone is not sufficient to induce the formation of a lattice, see Figs. 1b versus 1d. In system \scalerelS, a particle moves until it no longer feels its neighbours. At low concentrations, the system is quenched in a disordered state. Adding rotation to the particles, system \scalerelS, stirs and reshuffles them. Initially the amount of interactions increases, until a steady state configuration is reached in which particles are equally distanced. The whole lattice then rotates as a rigid body.

Increased rotor activity yields faster ordering. To measure order in the system we look at the structure factor, $S({\bf q})$, and at the two-dimensional bond-orientational order parameter, $\Psi_{6}^{j}=(1/n_{j})\sum_{i}e^{i6\theta_{ij}}$. $\Psi_{6}^{j}$ measures the orientation and degree of hexagonal order around particle $j$ Nelson (2002), the sum is over the nearest neighbours of particle $j$ as found from Delaunay triangulation Delaunay (1934), $n_{j}$ is the number of nearest neighbours, and $\theta_{ij}$ is the angle between the bond connecting particles $i$ and $j$ and the $x$ axis (which is an arbitrary reference). Let $\langle\cdot\rangle$ define averages over all ensemble particles. We define the average local and global order parameters to be $\langle|\Psi_{6}|\rangle$ and $|\langle\Psi_{6}\rangle|$, respectively. The global order parameter is close to zero for all initial random configurations, and is one for a perfect infinite hexagonal lattice. As shown in Fig. 2, purely repulsive interactions lead to slight increase in the order parameter, but do not result in an ordered array. Increasing the circulation, $\Gamma$, increases the order in the system, saturating the average local bond orientational parameter at a value of around 0.9 for large circulations. Repulsive interactions promote the formation of an ordered state but are unnecessary for maintaining order: once the system has reached the state shown in Fig. 1d, turning off the repulsive interactions does not destroy the order. This is a manifestation of the neutral stability of the purely hydrodynamic system, \scalerel*S.

Estimating the lattice rotation rate. Once formed, the crystals are not stationary in space, but rotate around their center of mass with an angular velocity $\Omega(r)$. Figure 3 shows the angular velocity as a function of radius for a finite system of particles interacting via Eq. 3 and with exponential repulsive interactions. We note two things about the angular velocity. First, it decreases with time (Fig. 3a). This decrease is a result of the collective expansion driven by the exponentially repulsive interactions. Second, changing the number of rotor proteins, but keeping the initial area fixed results in an increase in $\Omega$, as can be seen in Fig. 3b. In order to explain these features, we construct a simplified model based on two observations from Figures 1d and 3a: 1) the lattice rotates as a nearly rigid body, and 2) it forms a nearly perfect lattice (with spacing that only weakly depends on the radius). We hence assume rigid body rotation of a perfect lattice with spacing $d$ and angular velocity $\bar{\Omega}$. With these simplifications at hand, consider the following argument: for an infinite system the velocity must be zero from symmetry. To see this, consider the central particle in Fig. 3c (marked with a cross), each two opposing particles surrounding it will create opposite flows resulting in zero net velocity. Now note that for an infinite system all particles are identical, therefore the velocity of each particle in an infinite lattice is zero. For a finite system, it is thus the perimeter which dictates the angular velocity of the lattice, namely the number of proteins on the edge, $n_{r}$, and the distance between them, $d$, see Fig. 3c. Specifically, in the limit of small distances, $R\ll\lambda$, the velocity of the $k$th protein in complex notation, $z=x+iy$, is given by

[TABLE]

where distance is normalized by $\lambda$, and time by $2\pi\lambda^{2}/\Gamma$. Assuming the angular velocity is the same for all proteins, we can calculate the induced angular velocity of a single protein, for example, the particle just left of the center, $z_{k}=-De^{i\Omega t}$ (marked with a disk in Fig. 3c), where $D=d/\lambda$. The only components of the sum which do not cancel are the four rows at the edges (see Fig. 3). For a large number of particles, $n_{r}\rightarrow\infty$, the sum simplifies to,

[TABLE]

Similarly, in the opposite limit, of proteins separated by distances larger than the SD length, $\lambda$, we get $\Omega\propto 1/(D^{3}n_{r}^{3})$, (see supplementary information for more details). A more precise calculation, accounting for density and angular velocity variations as a function of radius, will lead to correction of the prefactor in Eq. 5, but the scaling will remain. After scaling the angular velocities at different times and for different system sizes according to Eq. 5, all results fall on a single curve, see Fig. 3d.

Effect of thermal noise. We have neglected thermal fluctuations so far but proteins are small ($\sim 10$nm) and prone to thermal forces. Thermal fluctuations must be added carefully when hydrodynamic interactions are included, as the thermal motion of one protein will effect that of another via the fluid in which they are immersed. To account for this, we use the positive definite analogue of the Rotne-Prager mobility tensor for a membrane, given in Ref. Sokolov and Diamant (2018). Figure 4 presents the results of increasing temperature on the global and local bond orientational order parameter for two different systems (of 200 and of 400 rotor proteins). At low temperatures the systems are ordered, transitioning into disordered ones at about $k_{\rm B}T/\tilde{\tau}\sim 0.3$, where $\tilde{\tau}={\eta_{2D}}\Gamma$ is the torque-dipole strength acting on a rotor protein.

To estimate the biologically relevant regime, we take the Péclet number to be ${\rm Pe}=va/D$, where $v$ is the hydrodynamic advective velocity of Eq. 3 summed over all particles in the ensemble, and $D$ is the diffusion coefficient of a rotor protein. Using measurements in the literature for membrane viscosity, protein rotation rate and size Kinosita et al. (1998); Levitt et al. (2009); Ueno et al. (2005), we estimate the Péclet number for ATP synthase to be around ${\rm Pe}\sim 2$. In our simulations, this corresponds to a normalized energy of about $k_{\rm B}T/\tilde{\tau}\sim 0.1$, which is below the transition point. This implies that biological systems should be in a relatively ordered state, but that thermal noise is also not negligible (see the middle snapshot in Fig. 4). Note however that this estimate does not incorporate the repulsive interactions, which would have their own Péclet number, and which we expect to further stabilize the system to thermal forces. Even in the absence of complete order, strong hydrodynamic interactions due to the rotation of the proteins may aid in mixing or in transporting other materials in the membrane.

Discussion. To conclude, a system of rotor proteins self organizes into a hexagonal lattice through a combination of hydrodynamic and repulsive interactions. At low rotor protein concentrations, with only repulsion, the system (\scalerelS) is quenched in a disordered state. Adding rotation (\scalerelS), and hence hydrodynamic interactions, shuffles the positions of the particles, driving them to an ordered steady state. In this sense rotation is similar to temperature, allowing the system to sample phase-space. Unlike temperature, the rotational interactions, in the absence of repulsive interactions, preserve the second moment as well as the Hamiltonian. The rotating system is then self-confined, whereas with temperature, an initial configuration will spread to infinity as the square root of time.

Throughout this work, we assumed the membrane to be completely flat. In the case of huge vesicles such as those used in Jiko et al. (2015), and for small thermal fluctuations, this assumption is reasonable. In a biological setting, however, ATP synthase is abundant in the highly curved mitochondrion. It was found that along the mitochondria ridges, ATP synthase dimers form a zigzag topology Blum et al. (2019), and that there is coupling between the membrane curvature and the activity of ATP synthase Almendro-Vedia et al. (2017). In the future, we plan to investigate the hydrodynamic coupling between rotation and curvature.

Acknowledgment. We thank Matan Yah Ben Zion, Enkeleida Lushi, and Florencio Balboa Usabiaga for useful discussions.