Dynamics of self-assembly of model viral capsids in the presence of a fluctuating membrane.

A coarse-grained computational model is used to investigate the effect of a fluctuating fluid membrane on the dynamics of patchy-particle assembly into virus capsid-like cores. Results from simulations for a broad range of parameters are presented, showing the effect of varying interaction strength, membrane stiffness, and membrane viscosity. Furthermore, the effect of hydrodynamic interactions is investigated. Attraction to a membrane may promote assembly, including for subunit interaction strengths for which it does not occur in the bulk, and may also decrease single-core assembly time. The membrane budding rate is strongly increased by hydrodynamic interactions. The membrane deformation rate is important in determining the finite-time yield. Higher rates may decrease the entropic penalty for assembly and help guide subunits toward each other but may also block partial cores from being completed. For increasing subunit interaction strength, three regimes with different effects of the membrane are identified.

Introduction

The formation of the protein shell of viruses has, due to its relative simplicity and importance in many diseases, become one of the most well-studied examples of self-assembly.[1] Although viruses are typically assembled within the cells of their host, the process may also be triggered in a bulk solution of viral proteins by changing the pH.[2] Such experiments have stimulated the application of simple computational models[1,3−7] to help understand assembly processes.

While much modeling has focused on the formation of virus capsids in the bulk, in recent work investigating the growth of viral shells around their genome, the assembly of simple subunits attracted to a flexible polymer was simulated.[8,9] Interaction with the polymer was found to allow assembly for parameters for which it would otherwise not occur. Encapsulation of spherical nanoparticles has also been considered both in experiment[10] and in simulation.[11,12] Experimentally, it was demonstrated that shells resembling different types of viral particles could be assembled by varying the nanoparticle diameter.

Beyond interactions with an encapsulated genome, there is also much evidence that membranes play an important role in assembly for many viruses.[13−21] In a recent publication,[22] we presented results on the effect of fluctuating membranes on the equilibrium of a system of self-assembling patchy colloids, designed to assemble viral core-like structures, from Monte Carlo (MC) simulations.[23] We found a nonmonotonic dependence of the promotion of assembly on membrane stiffness, as well as the formation of membrane buds. It is of course true that such effects would be observable in an analogous experimental system after sufficient time and to be expected that they will influence the products of dynamical assembly. However, on relevant time scales, self-assembly processes may not reach equilibrium and the products may be affected, for example, by kinetic traps.[1,24] It is therefore of foremost interest to consider simulations with realistic dynamics. Key dynamical features that we capture in our simulations are the viscosity of the membrane and hydrodynamic interactions, the inclusion of which may alter dynamics both quantitatively and qualitatively.[25]

Two key factors in the present work are attractions to the fluctuating membrane and hydrodynamic interactions. Previous computational studies have looked into the effects of each of these individually on the clusters formed by isotropic spherical colloids. Hydrodynamic interactions were found to change both the size and shape of clusters,[26] while attraction to a membrane was found to induce the formation of linear chains on the surface.[27] Further, attractions of particles to a membrane surface may cause the formation of buds[22,28,29] or tube-like structures.[30,31]

Here, as a simple model to gain insight into the effect of membranes on the dynamics of self-assembly, we consider primarily the same, patchy-particle, subunits,[6,32] which may assemble 12-component cores, as in our previous work,[22] and simulate their assembly using a dynamically realistic method. As previously, our subunits are coupled to a membrane modeled using particles bonded to form a triangulated surface.[33,34] The target core structure has icosahedral symmetry, similar to many viruses, although in reality enveloped viruses are larger. The remainder of the paper is organized as follows. In section 2 we describe our simulation models and in section 3 we present results from MC simulations on the equilibrium of the system. We then move on to dynamical simulations, describing simulation methods in section 4. We present results for the 12-component cores in section 5 and compare them to some results for other cores in section 6. Finally, we conclude in section 7.

Simulation Models

Rather than only considering enough subunits to form just one target structure as in our previous work,[22] we now simulate 180, allowing a maximum of 15 complete cores to be assembled. While it is expected that in experimental and biological situations it is also likely that a larger number of subunits will be available than required for one complete structure, this choice was additionally made for computational efficiency, so that, on a feasible time scale, although assembly of all possible cores may not occur, some complete cores will form. We simulate a membrane composed of 1156 particles. The simulation setup is sketched in Figure 1a.

Figure 1

(a) Simulation setup. Subunits, which are all identical, are rendered in yellow, with positions, but not extents, of patches for interactions with other subunits in red. Positions of patches for interactions with the membrane particles are in blue. The membrane is modeled as a triangulated surface of bonded particles. The particles forming the surface edge are confined to a frame region, which is located at a distance rframe from the periodic boundaries. In simulations with hydrodynamics, a stochastic rotation dynamics (SRD) solvent composed of point particles is included. Interactions between SRD particles are effected by first dividing the entire system into a grid of cells of side l0. (b) The radial part of the inter-subunit or subunit-membrane potential, U(r), with a well-depth ε is split into attractive (green) and repulsive parts (red). (c) The attractive part is multiplied by factors of the form Fang(θ), where θ is an angle that depends on the relative orientation of the interacting particles. (d) Sketch of momentum transfer between SRD particles in a cell: (i) Only particles within one cell interact. (ii) Velocities are subtracted from all particles such that the center of mass velocity is 0. (iii) All velocities are rotated, as signified by the heavy arrow, around a random axis, by a given angle. (iv) The subtracted velocities are added back on so that total momentum is conserved.

The interactions between subunits, ss, and between subunits and membrane particles, ms, are identical to those used in our previous work[22] but we describe the important features again here. The potentials are based on a Lennard-Jones form. As shown in eq 1, the potential is split into attractive, Uatt, and repulsive, Urep, parts. The interaction of two particles, i and j, separated by (i ≠ j), with orientations Ω and Ω, either both subunits or a subunit and a membrane particle, is given bywhere the forms of Uatt and Urep are shown in Figure 1b. γarea, γatt and γorient are dimensionless factors that take different forms for ss and ms interactions. For ss-interactions, γarea = γatt = 1 and, as depicted in Figure 1a, there are 5 patches on each subunit, which are arranged symmetrically around a single ms patch. The minimum of Uatt is set to −ε. γorient is used to control the patch width, and it has the form of a product of three functions of the form shown in Figure 1c, see also the Supporting Information. For the first two factors, the argument is the angle between the interacting patches and the center-to-center vector, . The parameters for determining patch width, see Figure 1c, are set to θ0 = θ1 = 0.2. In contrast, for the third factor, the argument is the angle between the projections of the membrane patches onto the plane perpendicular to and θ0 = θ1 = 0.4. The third factor represents the torsional stiffness of protein interactions.[6]

For ms-interactions, the minimum of Uatt is set to −εms. In these interactions, only the subunits are patchy, having one patch. Parameters for the one orientational function composing γorient, see Figure 1c, are θ0 = π/4 and θ1 = 0.2, and there is no penalty for subunits rotating around . Since, typically, assembling proteins will only be able to access one side of a membrane, we choose to make only one side of the membrane in our simulations attractive to subunits.[35] This is achieved by setting γatt = 1 if a subunit interacts with the “upper” side and γatt = 0 if it interacts with the “lower” side. The γarea factor is proportional to the area of surface that surrounds the interacting membrane particle. The length scale for ss-interactions is chosen as σss = 2.5l0 and the length scale for ms interactions is σms = 1.75l0, where l0 is a simulation unit of length. For the exact functional forms used in the ss and ms interactions, see the Supporting Information.

The membrane is modeled as in ref.[35] but we describe the key features again here. As depicted in Figure 1a, the membrane is composed of particles bonded to form a triangulated surface. To include membrane fluidity, MC moves that flip bonds between different particles are included.[34] The typical separation between bonded membrane particles is l0, maintained by a potential that has a flat central region but diverges at 0.67l0 and 1.33l0, see Supporting Information. We perform simulations in a box of size 45l0 × 45l0 × 45l0, giving a subunit number density within the range for which yield was found to weakly depend on concentration.[6] As in our previous work,[22] we consider a range of εss that, at equilibrium in the bulk, covers the crossover to complete assembly of all cores. Although approximately centered around the same εss value, for the larger number of cores considered here, the crossover is broader[36] and so a wider range of εss is used. The same range of εms as in ref (22) is considered, chosen to cover the crossover from freely diffusing to membrane-bound structures. The stiffness of our membrane is controlled by a parameter λb, through a potential, Ubend = λb(1 – ·), applied to all pairs of neighboring triangles in the surface, where and are the unit normal vectors of the triangles. We simulate using the three middle values from our previous work,[22] λb = √3kBT, 2√3kBT and 4√3kBT: at equilibrium, this covers the crossover from cores being able to cause budding of the membrane to them not being able to. As discussed in our previous work,[22] this range of bending stiffness is at the lower end of that expected for biological membranes. Given that in viral budding[37] intrinsic curvature is expected to be important, which is neglected in our model, the bending stiffness in our simulation is most relevant in terms of the cost of deformation.

Although our focus is on dynamical simulation, we first investigate the equilibrium of the system for comparison. For this purpose, we use MC simulations, employing a similar approach as in our previous work.[22] On the other hand, for molecular dynamics (MD) simulations, we include hydrodynamic interactions using a stochastic rotation dynamics (SRD) solvent,[38] a coarse-grained method in which the fluid is represented by point particles. SRD particle interactions are effected by dividing the system into a grid of cells, of side l0, at regular intervals and exchanging momentum by a rotation through a certain angle of velocities relative to the cell center of mass velocity. This procedure is shown schematically in Figure 1d. To understand the influence of hydrodynamic interactions, we also simulate using a method that neglects them, Langevin dynamics (LD), in which the effect of the solvent is represented by uncorrelated random, as well as drag, forces.[39,40]

To simulate a tensionless membrane, rather than box rescaling,[22] we use a new membrane boundary condition, recently introduced by us,[35] which is compatible with SRD. The edge of the membrane is attached to a square frame, with sides positioned at a distance rframe into the simulation box, as depicted in Figure 1a. For those triangles in the surface that have a side that forms part of the membrane edge, a bending potential of the same form as that between neighboring triangles is applied, except that the unit normal of the triangle is compared to a unit normal to the frame-plane. The distance rframe may increase and decrease during the simulation. To allow for deformation, the number of membrane particles bonded to the frame may also vary, with corresponding changes to the number of bonds in the bulk of the surface, Nb–bulk. For more details of the membrane boundary condition, see the appendix of ref (35), the functional form of the confining potential is also given in the Supporting Information. For consistency, this approach is also used in MC and LD simulations, in which, of course, the solvent is absent.

Results from Equilibrium Simulations

We first present, in Figure 2, results from MC simulations on the yield of complete cores, defined to be a cluster of 12 bonded subunits, each unit making 5 bonds to other cluster members. Two subunits are defined to be bonded if their interaction energy is < −0.25εss. Interaction strengths for different simulations lie on a grid from 0.12 to kBT in spacings of 0.08kBT for εms and from 4.5 to 8.02kBT in spacings of 0.32kBT for εss. Systems at different parameters were run in parallel using Multicanonical Parallel Tempering.[41] For each data point in Figure 2, approximately 4 × 109 attempted MC moves were performed, including about 4 × 104 Hybrid MC moves,[42] as well as Aggregate Volume Bias moves.[43] These were both found to significantly speed up relaxation. The largest error for a single data point was estimated to be about 0.6. Similarly to our results with one core,[22] for λb = √3kBT and λb = 2√3kBT at higher εms, the assembly of the cores causes the membrane to form buds, although these now generally contain multiple cores. For λb = 4√3kBT, budding did not occur. Again as for single cores, for high εms, assembly occurs for lower values of εss: the membrane promotes assembly. Here, membrane-dependent, low εss assembly does not occur to the same extent as for high εss because, due to steric repulsion, only a fraction of the cores may interact with the membrane at once, typically about 4 cores in the case where a bud is formed. Whereas for one core, the range over which promotion occurred was clearly largest for λb = √3kBT, here the results for λb = 2√3kBT are very similar. This may be because multiple cores together effectively form a larger object deforming the membrane.

Figure 2

Results from MC simulations. Average yield of complete cores, ⟨yield⟩, as a function of subunit-membrane interaction strength, εms, and inter-subunit interaction strength, εss, for different membrane stiffnesses, λb: (a) λb = √3kBT; (b) λb = 2√3kBT; (c) λb = 4√3kBT.

Dynamical Simulation Methods

We next give details of our dynamical simulation methods. The SRD particles have mass m and number density per cell γ = 5. We define our unit of time, t0 = l0(m/kBT)1/2. Collisions are performed every Δtcoll = 10–1t0 and we use an SRD rotation angle of π/2, giving a fluid viscosity of ηf = 2.5m/l0t0.[44] We apply a SRD-cell level thermostat that conserves momentum to maintain the temperature.[38] Membrane particles are coupled to the SRD solvent by including them in the collision step.[38] There will typically be about one membrane particle per SRD cell and we set their mass to γm, giving a short-time friction coefficient ζmem = 15.8(m/t0).[44]

Unlike membrane particles, subunits have rotational degrees of freedom and so are coupled to the SRD solvent using bounce-back boundary conditions.[45] For their interactions with the fluid, subunits are treated as solid spherical particles of radius a = l0, having mass M = (4/3)πa3mγ and moment of inertia I = (2/5)Ma2. Every Δtbound = 10–2t0, the SRD particles are checked. If an overlap with a subunit is detected, then the SRD particle with velocity is first moved by (−1/2)Δtbound and then shifted radially to the edge of the subunit, from the center, where || = a. This scheme is based on the fact that for SRD particles the average crossing of the subunit boundary is halfway through a time step. It was found to function well in previous work.[46] A bounce-back collision is then performed: the radial, ⊥, and tangential, ∥, components of are updated according toHere, A = 2M/(m+M), B = 7m/2M and the surface velocity = + × , where is the center of mass velocity of the subunit and is its angular velocity around an axis that passes through the center of mass. Equation 2 is valid for I = 2/5Ma2. After all overlapping SRD particles have been rebounded, corresponding changes to the subunit velocity and angular velocity, Δ = (m/M)Σi(old – new) and Δ= (m/I)Σ × (old – new), where i indexes the different rebounded particles, are applied so that momentum and energy are conserved. If M ≫ m, SRD particle velocities relative to the surface are completely reversed; for our parameters, M ≈ 20m.

Overlapping of embedded particles in an SRD fluid may lead to a spurious depletion attraction.[47] In fact, even if particles are prevented from overlapping, the bounce-back scheme may need to be iterated due to the possibility of a fluid particle interacting with more than one solute particle within Δtbound. We avoid these issues by choosing the excluded volume length for subunit interactions, σss = 2.5l0, so that the typical closest approach of two subunit fluid surfaces ≈ 0.5l0 is much greater than the typical displacement of a fluid particle ≈ 10–2l0.

Bounce-back interactions between SRD particles and embedded colloids lead to spurious slips at the colloid surface. Methods exist to ameliorate this by the introduction of virtual particles but, for mobile colloids, this was found to lead to deviations from expected thermal distributions.[45] In our simulations, the concern is moot anyway, because of the discrepancy between the radii chosen for inter-subunit and subunit-fluid interactions. Effectively, there is a slip-velocity at the subunit surface, as defined by its interactions with other subunits, which has contributions from these two different sources. Given that the subunits are typically representing protein complexes, which are rough on length-scales up to many solvent molecules,[48] rather than smooth colloids, this is reasonable.

For bounce-back boundaries, the short-time friction coefficients for the subunits may be calculated using a modified Enskog theory.[45] For our parameter choice, this gives coefficients of ζ = 62.0(m/t0) and ζω = 73.3(ml02/t0), for linear and angular velocities respectively. Comparing the corresponding correlation times, M/ζ and I/ζω, to typical thermal velocities, we obtain values of 0.07l0 and 0.04 for the typical length and angular displacements over which the subunit motion is correlated. These are smaller than the typical separation of subunits and patches, given by σss = 2.5l0 and ≈ 0.4 respectively, so that at the scale that assembly occurs on, subunit motion is diffusive. Similarly, the length scale over which membrane particle motion is correlated is 0.1l0.

To obtain parameters for simulations without hydrodynamic interactions, we simulated single subunits and membrane particles in a box of the same size as that used for assembly with an SRD solvent. The friction coefficients extracted were lower than the short-time values due to long-time hydrodynamic contributions. These friction coefficients were input to LD simulations. In this way, the hydrodynamic contribution to the self-diffusion coefficients is included but hydrodynamic interactions between different particles are neglected. An alternative approach to simulating without hydrodynamics is to use an SRD fluid and randomize particle velocities at every step. For colloids, however, this has been found to introduce an unphysical caging effect.[49]

Membrane fluidity is included by performing a certain number of attempts to flip bonds between neighboring pairs of membrane particles[34] every 10–1t0. The membrane viscosity is set by the level of attempted bond-flips and we consider three different rates: Nb–bulk, 10–1Nb–bulk and 10–2Nb–bulk attempted bond-flips per 10–1t0, where Nb–bulk is the number of bonds in the bulk of the membrane and the resulting numbers are rounded to integers. By considering Poiseuille flows in two-dimensional membranes,[34] the corresponding membrane viscosity, ηm, may be estimated. For the highest rate of flips the value is estimated to be 35.1 ± 0.1m/t0,[35] whereas for the lower rates we estimate 133.3 ± 0.6m/t0 and 1190 ± 60m/t0 respectively. For a lipid bilayer in water, the ratio of membrane to fluid viscosities, lη, is typically around 1–10 μm.[50] In our simulations the solvent viscosity ηf = 2.5m/l0t0 so that, if our subunits represent capsomers with a size on the order of 10 nm,[51] then the ratio of their hydrodynamic radius to lη is around the expected range.

Results from Dynamical Simulations

We next present results from our dynamical simulations. Interaction strengths for dynamical simulations were chosen to coincide with those for MC simulations, although fewer were considered due to higher computational costs. A closer spacing between the highest interaction strengths was chosen as it was expected that the most interesting results would be found here. All averages are taken over at least 5 independent runs and in some cases over 10. We consider the same values of λb, and also εm and ε in the same range, as for the equilibrium MC simulations. We simulate primarily using SRD but, for λb = √3kBT and 2√3kBT, we also simulate using LD for comparison, to gain insight into the importance of hydrodynamic interactions. LD simulations were essentially identical to the SRD ones expect that, rather than having regular interactions with an explicit fluid, subunits and membrane particles were, at each MD integration step, subject to random and friction forces.[39,40] The system was initially simulated for either 8 × 103t0 with SRD, or 2 × 104t0 with LD, without attractive interactions. These times were chosen as being sufficient to allow membrane relaxation. Subsequently, attractions were switched on and the system was simulated for a further 5 × 104t0 to gather results. In contrast to the MC results, for the stiffest membrane, λb = 4√3kBT, with the highest membrane-subunit interaction strength only, εms = kBT, in some, though not all runs, budding occurred.

Were it possible to run the dynamical simulations indefinitely, it is expected that results would eventually converge to those found for the equilibrium simulations. However, as the simulation progresses, further assembly becomes increasing slow as the supply of free subunits is depleted and eventually relies on rearrangement of subunits between partially formed structures, possibly moving into or out of a membrane bud. It is thus necessary to choose a finite simulation time shorter than that required for complete assembly and inevitably the results obtained will depend on it. For our chosen simulation time, the highest yield observed in any simulation is ∼50% of the possible maximum. It is nonetheless sufficient for the effect of the membrane on assembly to be apparent. However, the finite time chosen should be borne in mind when considering the results.

First, in Figure 3, we plot the averages of various quantities as a function of ε and ε at different simulation times, t, for λ = 2√3kBT and ηm = 133.3m/t0. Considering Figure 3a, at later times, the largest number of correctly assembled cores are obtained for the second highest εss, 7.38kBT. This is close to the optimal value obtained in previous work[6] with a very similar model of about 7.14kBT. Although for the highest value, εss = 8.02kBT, the total interaction energy between subunits relative to the interaction strength is somewhat lower, see Figure 3b, this corresponds to many incomplete cores assembling, thus starving the system of free subunits. This kinetic trap is not related to the membrane and has often been observed previously.[1] In contrast, for high εss, increasing attraction to the membrane hinders complete assembly due to the membrane enveloping, or partly surrounding, partial cores too quickly, thus preventing subunits or other partial cores from approaching them. The fast envelopment is apparent in Figure 3, where it may be seen that, for high ε and ε, the interaction energy of the subunits with the membrane approaches its final value much more quickly than the yield. Similarly to the MC results, a promotion of the finite-time assembly for low εss at high εms occurs.

Figure 3

Plots as a function of subunit-membrane interaction strength, εms, and inter-subunit interaction strength, εss at different times, t, increasing from left to right in intervals of 1 × 104t0, for membrane stiffness, λb = 2√3kBT and membrane viscosity ηm = 133.3m/t0, from SRD simulations. (a) Average yield of complete cores, ⟨yield⟩. (b) Average total interaction energy between subunits, relative to the interaction strength, ⟨Uss/εss⟩. (c) Average total interaction energy between subunits and the membrane, relative to the interaction strength, ⟨Ums/εms⟩.

Comparing results for equilibrium, Figure 2, to those from dynamical simulations, there are several differences. Clearly, for equilibrium results, kinetic traps do not play a role. Furthermore, the interaction strength at which assembly starts is slightly lower at equilibrium than after a finite time in dynamical simulation. For the lowest εms, εms = 0.12kBT, where the membrane does not play a significant role and assembly occurs in the bulk, whereas at equilibrium there are complete cores at εss = 5.46kBT, albeit at a relatively small yield, no complete cores were formed within the allowed time in dynamical simulations. Similarly, the range of parameters for which there is assembly promotion is larger for the equilibrium results.

Considering Figure 3a, we note that a time lag before complete cores are assembled, seen in previous work,[1] is apparent for many data points at t = 104t0, including for εss = 7.38kBT and low εms, where the yield is highest at later times. However, for some data points, primarily with high εms, some complete cores are already present at t = 104t0: as well as causing a higher yield for low εss once assembly has progressed significantly, attraction to the membrane may also speed up the formation of a single core. By confining subunits to a surface, the effective size of the space that they must search to find each other is reduced. The membrane may also mediate effective attractions, directing subunits and partial cores toward each other[28] and, if deformation occurs, it may bring membrane-attached subunits closer together. Conversely, deformation of the membrane may also tend to block assembly, preventing partial cores from being accessed by subunits or other partial cores, leading to an increase in assembly time. The effect of the membrane on single-core assembly times is also shown in Figure 4, where we plot the average time until the first complete core in the system is formed, which we denote ⟨τ⟩, as a function of εms. We note that, since this quantity is based on a single assembly event, large fluctuations were seen for lower interaction strengths and for some parameters additional simulations were run. For both membrane stiffnesses shown, λb = 2√3kBT and 4√3kBT, ⟨τ⟩ tends to be lower for high εms for εss = 6.42kBT, whereas for εss = 7.38kBT the curve is flatter. For subunit interaction strengths that are approximately optimal for bulk assembly, the process is sufficiently fast that the membrane does not affect ⟨τ⟩, whereas for lower values it may cause a significant speed-up. However, we note that, for high εms with some parameters, ⟨τ⟩ shows an increase. This is consistent with the membrane blocking assembly.

Figure 4

Average time until the first complete core is assembled, ⟨τ⟩, as a function of membrane-subunit interaction strength, εms, for inter-subunit interaction strength, εss = 6.42kBT (red) and εss = 7.38kBT (green), with different membrane viscosities: ηm = 35.1m/t0 (▲); ηm = 133.3m/t0 (●); ηm = 1190m/t0 (■). (a) Membrane stiffness, λb = 2√3kBT. (b) λb = 4√3kBT.

We next consider the average yield of complete cores, ⟨yield⟩, measured at the end of the simulation, at time t = 5 × 104t0. We present results for all the different parameter sets we have simulated, except for the lowest εss, for which only a very small amount of assembly occurred at the highest εms. We split the results into 3 figures by membrane stiffness: λb = √3kBT in Figure 5; λb = 2√3kBT in Figure 6; λb = 4√3kBT in Figure 7. Within each figure, results are divided by εss into four different subfigures labeled a–d. ⟨yield⟩ is plotted as a function of εms, with curves corresponding to different membrane viscosities and simulation methods indicated by different colors, symbols and line types. Rather than describe in detail the specific features of each figure, we discuss the general trends that arise and pick out interesting features.

Figure 5

Plots of the average yield of complete cores, ⟨yield⟩, as a function of membrane-subunit interaction strength, εms, for membrane stiffness, λb = √3kBT at t = 5 × 104t0. From simulations with SRD (solid lines, filled symbols) or LD (dashed lines, open symbols), for different membrane viscosities: ηm = 35.1m/t0 (blue, ▲/Δ); ηm = 133.3m/t0 (green, ●/○); ηm = 1190m/t0 (red, ■/□); and subunit interaction strengths, εss, as indicated in panels a–d.

Figure 6

Plots of the average yield of complete cores, ⟨yield⟩, as a function of membrane-subunit interaction strength, εms, for membrane stiffness, λb = 2√3kBT at t = 5 × 104t0. From simulations with SRD (solid lines, filled symbols) or LD (dashed lines, open symbols), for different membrane viscosities: ηm = 35.1m/t0 (blue, ▲/Δ); ηm = 133.3m/t0 (green, ●/○); ηm = 1190m/t0 (red, ■/□); and subunit interaction strengths, εss, as indicated in panels a–d.

Figure 7

Plots of the average yield of complete cores, ⟨yield⟩, as a function of membrane-subunit interaction strength, εms, for membrane stiffness, λb = 4√3kBT at t = 5 × 104t0. From simulations with SRD, for different membrane viscosities: ηm = 35.1m/t0 (blue, ▲/Δ); ηm = 1190m/t0 (red, ■/□); and subunit interaction strengths, εss, as indicated in panels a–d.

Overall, as εss is increased, we identify three different trends for finite-time assembly with realistic dynamics. First, at εss = 5.46kBT and 6.42kBT, both lower than the bulk-assembly optimal value, increasing εms, which tends to increase the membrane deformation rate in all regimes, may promote finite-time assembly. With high εms, assembly also occurs for values of εss where there is no bulk assembly within our simulation time. However, the rate at which the membrane deforms is important. It is influenced by various factors, for example the strength of the attraction of subunits to the membrane or membrane viscosity, and results in competing effects on assembly. Increasing it via εms may at first aid assembly, as seen in the initial increase in ⟨yield⟩ with εms in parts a and b of Figures 5–7, but when it is too high, the yield may decrease again, as is seen particularly clearly in Figure 6a and b. Results depend on factors such as membrane viscosity and hydrodynamic interactions: decreasing membrane viscosity or including hydrodynamic interactions both increase the deformation rate. Attraction to, deformation of, and encapsulation within, the membrane is expected to decrease the subunit entropy such that the difference in entropy between unassembled and assembled states is less. Additionally, the deformation of the membrane may help to guide subunits attached to it toward each other. It is expected that these effects will all play a role in promoting assembly, although their relative importance may not be easily deduced. However, if the deformation occurs too quickly, before complete cores are formed, the membrane will hinder further subunits, or other partial cores, from approaching the partial structure, preventing its completion.

Interestingly, in this low ε regime, finite-time assembly is promoted even for λb = 4√3kBT, Figure 7a,b, although for this membrane stiffness budding only occurs for the highest εms. For this membrane stiffness, results do not depend on membrane viscosity, confirming that here budding does not play a role. Despite a lack of envelopment, it is expected that attachment to the membrane will nonetheless reduce the entropic cost of forming a partial core. A second plausible mechanism is a local increase in subunit density near the membrane surface. In contrast, for the lower two λb, Figure 5a,b and Figure 6a,b, the membrane deformation rate does play a role. Comparing SRD and LD results, simulations with hydrodynamic interactions show a larger promotion of finite-time assembly as εms is increased, at least initially. For some parameters the yield decreases again as εms is increased further and this drop off occurs earlier with hydrodynamic interactions. Furthermore, particularly for λb = 2√3kBT, membrane viscosity, ηm, is also important. Especially for SRD results, decreasing η shifts the point at which the finite-time yield begins to decrease to lower εms. Interestingly, the effect of hydrodynamic interactions and membrane viscosity are much stronger for λb = 2√3kBT than for λb = √3kBT.

In the second regime, for εss = 7.38kBT, at about the bulk-assembly optimal value, increasing εms tends to decrease the finite-time yield. Here, there is no strong dependence on membrane viscosity or hydrodynamic interactions and, additionally, results are quite similar for all three membrane stiffnesses. This suggests membrane deformation is not crucial, rather the decrease in yield may occur because attraction to the membrane promotes the faster assembly of partial cores, bringing the system into the monomer starvation trap that is only seen in the bulk for higher εss.

Finally, for the highest εss, εss = 8.02kBT, where there is a monomer starvation kinetic trap for bulk assembly, there is no clear effect of the membrane on finite-time assembly. Since assembly, at least of partial cores, occurs very quickly in the bulk, results here are likely dominated by nonmembrane-associated assembly.

In Figure 8, we show snapshots of final configurations from simulations with λb = 2√3kBT and ηm = 133.3m/t0. For εms = kBT and εss = 7.38kBT, Figure 8a, for which the average yield was reduced compared to the low εms value, three of the four cores encapsulated in a bud that are depicted are incomplete. This snapshot corresponds to the optimal εss regime. In contrast, for εms = 0.6kBT and εss = 6.42kBT, Figure 8b, for which the average yield was enhanced compared to low εms value, only one of the four cores encapsulated, or partially encapsulated, in a bud that are depicted is incomplete. This snapshot corresponds to the low εss regime. Although these snapshots only depict the situation in two individual runs, they illustrate how the membrane may block assembly completion when its deformation rate is too high by preventing partial cores from being accessed by subunits or other partial cores.

Figure 8

Snapshots from simulations with SRD at t = 5 × 104t0 with membrane stiffness, λb = 2√3kBT and membrane viscosity, ηm = 133.3m/t0: (a) Membrane-subunit interaction strength, εms = kBT and inter-subunit interaction strength, εss = 7.38kBT; (b) εms = 0.6kBT, εss = 6.42kBT. Subunits are shown in yellow and membrane particles in green. Only subunits within 6l0 of a membrane particle are plotted. Membrane particle size has been reduced to make structures within buds more visible. Completed cores are circled in blue, while partially assembled ones are circled in red.

A useful quantity to indicate the amount of membrane deformation is rframe, the distance from the edge of the simulation box of the frame to which the edge of the membrane is bound, which increases as the membrane distorts its shape out of the plane. To show how changing membrane viscosity, and including hydrodynamic interactions, alters the rate and extent of membrane deformation, we plot, in Figure 9, ⟨rframe⟩ as a function of time with λb = 2√3kBT and εss = 6.42kBT for different εms from 0.36kBT. Apart from the lowest εms, the membrane deformation occurs faster and to a greater extent for simulations with hydrodynamics. Hydrodynamic interactions increase the rate of budding. Since budding requires the whole of the membrane to move, correlations mediated by hydrodynamics promote it. Furthermore, for SRD simulations at the highest two εms, there are also significant differences between membrane viscosities with a trend as expected: ⟨rframe⟩ is largest for the smallest viscosity.

Figure 9

Plots of the average position of the frame, ⟨rframe⟩, as a function of time, t, for membrane stiffness, λb = 2√3kBT and inter-subunit interaction strength, εss = 6.42kBT. From simulations with SRD (solid lines, filled symbols) or LD (dashed lines, open symbols), for different membrane viscosities: ηm = 35.1m/t0 (blue, ▲/Δ); ηm = 133.3m/t0 (green, ●/○); ηm = 1190m/t0 (red, ■/□); and subunit-membrane interaction strengths: (a) εms = 0.36kBT; (b) εms = 0.6kBT; (c) εms = 0.84kBT; (d) εms = kBT.

We next consider the distributions of cluster size, nc, in our SRD simulations, H(nc). In Figure 10, we plot distributions for λb = 2√3kBT and ηm = 35.1m/t0 at three different times. For εss = 6.42kBT, the lowest εss for which, in dynamical simulations, there is bulk assembly, there are few clusters with intermediate sizes when εms is low. At all times considered, the majority of clusters are of size two; at later times there is an additional peak at size 12, corresponding to complete cores. In contrast, when εms is high, the attraction to the membrane stabilizes intermediate cluster sizes at early times. At later times, the clusters have grown but the peak near 12 is less sharp, with similar numbers of cores of size 10 and 11, and also some larger ones. This shows the effect of the membrane blocking the completion of partial cores. It also seen for higher εss, εss = 7.38kBT, where the distribution for low εms is much flatter at early times with many clusters of intermediate sizes.

Figure 10

Results from SRD simulations. Histograms, H, of cluster size, nc, for membrane stiffness, λb = 2√3kBT and membrane viscosity, ηm = 35.1m/t0 at different times, t: t = 1 × 104t0 (red); t = 3 × 104t0 (green); t = 5 × 104t0 (blue); with different parameters: (a) Membrane-subunit interaction strength, εms = 0.12kBT and inter-subunit interaction strength, εss = 6.42kBT; (b) εms = kBT, εss = 6.42kT; (c) εms = 0.12kBT, εss = 7.38kBT; (d) εms = kBT, εss = 7.38kBT.

In Figure 11, we show the effect of membrane viscosity on the cluster size distribution. We plot results from the end of the simulations, at time t = 5 × 104t0, again for λb = 2√3kBT. With εss = 6.42kBT, we see that, for the highest viscosity, increasing εms leads to a distribution that is more strongly peaked at 12. In contrast, for the lowest viscosity, although increasing εms does lead to more larger clusters, it also gives a much broader distribution around 12. A similar effect occurs for εss = 5.46kBT: for this εss too, when membrane viscosity is low, high εms causes the membrane to encapsulate the assembling capsids too quickly, blocking their completion.

Figure 11

Results from SRD simulations. Histograms, H, of cluster size, nc, for membrane stiffness, λb = 2√3kBT at time, t = 5 × 104t0 for different membrane-subunit interaction strength values: εms = 0.12kBT (red); εms = 0.6kBT (green); εms = kBT (blue); with different parameters: (a) Membrane viscosity, ηm = 1190m/t0 and inter-subunit interaction strength, εss = 6.42kBT; (b) ηm = 35.1m/t0, εss = 6.42kBT; (c) ηm = 1190m/t0, εss = 5.46kBT; (d) ηm = 35.1m/t0, εss = 5.46kBT.

As in similar previous models,[5] the concentration of our subunits is relatively high compared to experimental systems, and furthermore the number of subunits in a completed core is low. These choices are necessary for computational tractability but have the consequence that the assembly rates in our simulations are much higher than experimental ones. Assuming that subunits correspond to capsomers of size on the order of 10 nm, and matching the drag coefficient of our subunits, we estimate that our simulation length is around 5 ms, whereas in vitro(24) and in vivo(52) experiments have observation times on the order of minutes. Thus, a direct quantitative comparison cannot be made. Our results rather demonstrate how the rate of membrane deformation compared to the assembly rate may affect the success of the latter. Furthermore, they show how properties such as membrane viscosity, which might be varied experimentally by changing lipid composition[53] are expected to impact on the assembly process. The interactions between assembling viral capsomers are typically relatively weak[54] as this aids the avoidance of kinetic traps and thus the first, low ε regime identified is most likely to be relevant to these systems.

Other Target Cores

We have investigated the effect of a membrane on core assembly of icosahedral cores. It is expected that many of the qualitative features of the results, such as the interplay between the membrane promoting assembly by confining subunits and hindering it by blocking additional subunits or other partial cores from approaching partial structures, will be general to other target structures. To gain more insight into the transferability of the findings to other core shapes, we finally, in Figure 12, present results on the yield obtained in SRD simulations when the target structure is changed from the icosahedral core. As for icosahedra, we define a cluster to be a complete structure when it contains the correct number of subunits for the target structure and each subunit patch on each subunit forms a bond with another member of the cluster. We choose subunits with patches such that their interactions are minimized for cubic and dodecahedral structures.[6] Otherwise, parameters such as patch width are unchanged, with the membrane patch still lying on the symmetry axis as defined by the subunit patches and pointing outward in a complete structure. In both cases, the subunits only have three patches for bonding with other subunits and thus form fewer bonds in a complete structure than the icosahedral subunits. Correspondingly, the range of εss was shifted up by about 2kBT but the range of εms remained the same.

Figure 12

Results from SRD simulations. Plots of the average yield of complete structures, ⟨yield⟩, as a function of membrane-subunit interaction strength, εms and inter-subunit interaction strength, εss at time t = 5 × 104t0 for membrane stiffness, λb = 2√3kBT and membrane viscosity, ηm = 133.3m/t0: (a) Subunits with interactions to form a cube. (b) Subunits with interactions to form a dodecahedron. Note the different scales for ⟨yield⟩ and also the higher values of εss as compared to the results for icosahedral cores.

We simulated for λb = 2√3kBT and ηm = 133.3m/t0. For both cubes and dodecahedra, as for icosahedra, budding of the membrane occurred for high εms. In Figure 12a, many of the features of the finite-time assembly of icosahedra are reproduced for cubes. For high εms, there is finite-time assembly at lower εss than for low εms. There is also an increase in finite-time yield with increasing εms for the lowest εss for which assembly occurs without significant attraction to the membrane. Unlike for icosahedra, at least for these membrane parameters, the yield does not drop off again as εms is increased further. This may be because, since they are composed of less subunits, cubes assemble faster than icosahedra. For the highest εss, however, there is a reduction in finite-time yield with εms, similar to results for icosahedra.

As observed in previous work,[6] finite-time yields of dodecahedra, Figure 12b, were low. However, here again, there is evidence that attraction to the membrane may promote assembly for εss for which it would otherwise not occur. Although it is not apparent in our results, as seen in previous work with a very similar model,[6] it is expected that, if the subunit interaction strength were increased sufficiently, the same nonmembrane-related kinetic trap that is observed for icosahedra would also be seen for cubes and dodecahedra.

Conclusions

To summarize, we have applied a simple patchy-particle model to investigate the effect of interactions with a fluctuating membrane on the dynamics of the assembly of core structures with the same symmetry as many viral cores. As well as interaction strengths, the key parameters we varied were membrane stiffness and viscosity. We also considered the effect of hydrodynamic interactions by simulating both with SRD and LD. As at equilibrium, for assembly with realistic dynamics, attraction to a membrane may promote finite-time assembly, also for subunit interaction strengths, εss, for which it does not occur in the bulk. Furthermore, for ε less than the optimal bulk value, attraction to the membrane also decreases the single-core assembly time.

Membrane budding occurred in dynamically realistic simulations and its rate was strongly increased by hydrodynamic interactions, as well as by lowering the membrane viscosity. The rate of membrane deformation is important in determining the assembly yield after finite time. Relatively high rates may promote assembly by increasing the envelopment of assembling cores and thus decreasing the entropic penalty and also by guiding subunits toward each other. However, if the rate is too high, the membrane may block partial cores from being completed. Three regimes with different effects of the membrane were identified. For εss less than the bulk optimum, finite-time yields depend intricately on a combination of all parameters and may both increase and decrease as attraction to the membrane is increased. For εss about equal to the bulk optimum, finite-time yields do not depend strongly on the membrane deformation rate and tend to decrease as attraction to the membrane is increased. For εss higher than the bulk optimum, assembly in the bulk is affected by a monomer starvation kinetic trap and the membrane has little influence.

Finally, results with qualitative similarities were also found for core structures with cubic and dodecahedral symmetries. In future work it would be interesting to investigate more different structures, in particular much larger cores.