## Abstract

Retrovirus particle (virion) infectivity requires diffusion and clustering of multiple transmembrane envelope proteins (Env_{3}) on the virion exterior, yet is triggered by protease-dependent degradation of a partially occluding, membrane-bound Gag polyprotein lattice on the virion interior. The physical mechanism underlying such coupling is unclear and only indirectly accessible via experiment. Modelling stands to provide insight but the required spatio-temporal range far exceeds current accessibility by all-atom or even coarse-grained molecular dynamics simulations. Nor do such approaches account for chemical reactions, while conversely, reaction kinetics approaches handle neither diffusion nor clustering. Here, a recently developed multiscale approach is considered that applies an ultra-coarse-graining scheme to treat entire proteins at near-single particle resolution, but which also couples chemical reactions with diffusion and interactions. A model is developed of Env_{3} molecules embedded in a truncated Gag lattice composed of membrane-bound matrix proteins linked to capsid subunits, with freely diffusing protease molecules. Simulations suggest that in the presence of Gag but in the absence of lateral lattice-forming interactions, Env_{3} diffuses comparably to Gag-absent Env_{3}. Initial immobility of Env_{3} is conferred through lateral caging by matrix trimers vertically coupled to the underlying hexameric capsid layer. Gag cleavage by protease vertically decouples the matrix and capsid layers, induces both matrix and Env_{3} diffusion, and permits Env_{3} clustering. Spreading across the entire membrane surface reduces crowding, in turn, enhancing the effect and promoting infectivity.

This article is part of the themed issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

## 1. Introduction

Retrovirus particles (virions), in particular of HIV-1, are complex, crowded, macromolecular entities formed by an array of striking morphological changes [1,2]. Viral polyproteins, Gag and GagPol, containing sequentially connected matrix (MA) then capsid (CA) at their N-termini (figure 1), assemble at the intracellular host membrane as a tightly packed hexameric-holed lattice [3,4] mediated by both intra- and inter-hexamer CA interactions [5]. The lattice is monotopically anchored to the membrane by the myristoylated N-termini of MA, where MA forms trimers [6] and has been observed *in vitro* to form hexamers of trimers [7] in planar assemblies. Budding leads to nascent spheroidal virions of varying size (radius ∼63±5 nm) and results in complete envelopment by, but incomplete lattice coverage (approx. 64±10%) of, host membrane [8]. Subsequent maturation involves autocatalysed chemical degradation of Gag & GagPol [9] by viral protease (PR) [10] coupled to self-assembly of liberated CA [11], leading to a conical CA shell [12–14] that encapsidates ribonucleoprotein and viral enzymes (figure 1) [15].

The CA assembly pathway remains incompletely understood. One view holds that liberated CA condense to form the mature capsid [16,17] through a displacive pathway that maintains the lattice intact but remodels its shape. This has been recently supported by observations of non-diffusional rolling of an intact CA lattice [18] and unclosed curled sheets during assembly [19]. Another view holds that capsid subunits disassemble upon cleavage and then reassemble to form the conical capsid in a two-step process [20,21], based on observations (i) that the immature CA lattice is more dense [20], (ii) that a significant fraction of virions form two capsid cores [22], and (iii) that mature virions are composed of single shells with about 1500 CA molecules (substantially fewer than the total number (approx. 2400) of Gag [8]) and a significant pool of unassembled CA subunits [22,23].

Infectious virions are capable of membrane fusion with host cells—mediated by the viral transmembrane envelope glycoprotein trimer (Env_{3}) [24]. Env_{3} is composed of a major ectodomain (ED), a transmembrane domain (TMD) and a cytoplasmic tail (CT) [25] (figure 1). The latter extends into the MA layer, interacts with both the membrane and MA and is most likely situated within the hexameric holes of the Gag lattice [24]. Typical virions contain 5–15, initially separated Env_{3} molecules [26]. Productive infection often requires and is enhanced by oligomerization of multiple Env_{3} [27]. Activation of fusion activity is triggered by PR-dependent Gag lattice degradation [28] and concomitant with both internal virion remodelling and surface diffusivity of Env_{3} [29].

Spatio-temporal *in vivo* characterization of the coupling between Gag degradation, capsid assembly and Env_{3} oligomerization is difficult to examine in macromolecular detail. Computational simulations can, in principle, inform such spatio-temporal resolutions. Previous approaches have focused on either (i) capsid self-assembly using both equilibrium [30–32] and non-equilibrium [33,34] coarse-grained molecular dynamics (CGMD), Monte Carlo [35] approaches and continuum models [19,36,37] or (ii) polyprotein degradation using reaction kinetics (RK) models [38,39]. All-atom MD has also informed large-scale virion capsid dynamics [14] but cannot access assembly time scales. Modelling virus assembly, including both icosahedral and conical HIV-1 capsid self-assembly, has been reviewed recently [40,41].

The spatio-temporal linkage between internal virion Gag degradation and diffusional oligomerization of Env_{3} has not been modelled and remains unclear. To do so requires development of a coupled model of the mesoscopic-scale chemical reaction, diffusion and oligomerization processes involved. CGMD and Brownian dynamics (BD) approaches handle diffusion and association [42] but exclude chemical reactions [43], thus the relation between internal remodelling due to polyprotein cleavage and Env_{3} diffusion is inaccessible. RK approaches can account for such cleavage events but not for diffusion, possible crowding phenomena or stochastic fluctuations due to small copy numbers of key particles.

A novel detailed modelling approach has been developed, termed interacting particle-based reaction–diffusion (iPRD) [44]. This combines space-excluded particle-based isotropic BD diffusion with state-changing chemical reactions, including the assignment of inter-particle and particle-geometry interaction potentials. Particle diffusion is carried out therein by numerical integration of the isotropic BD equation:
1.1where **x**(*t*) is the instantaneous position of a particle at time *t*, ∇*V* (**x**(*t*)) is the gradient of an applied interaction potential, *k*_{B} is the Boltzmann constant, *T* is the temperature, *η*(*t*) is a three-dimensional Wiener process (an independent random process with normally distributed increments) and *D* is the diffusion constant.

State-changing bimolecular reactions are based on a partition of the second-order macroscopic reaction rate constant (*k*^{mac}) into a Smoluchowski encounter rate constant (*k*^{enc}=4*πD*_{AB}*R*_{I,AB}, where *D*_{AB} = *D*_{A} + *D*_{B}, *D*_{A} and *D*_{B} are the diffusion constants of two reacting particles *A* and *B* and *R*_{I,AB} is the sum of the individual interaction radii, *R*_{I,AB} = *R*_{I,A} + *R*_{I,B}, that corresponds to an encounter) and a microscopic first-order rate constant (*k*^{mic}), obtained by fitting known or given macroscopic constants to the Erban–Chapman equation [45]:
1.2

Ultra-coarse-graining (UCG) enables entire proteins to be described at single or near-single particle resolution with simulation time steps on the nanosecond scale [46]. This spatio-temporal resolution is compatible with that of modelling lateral multiprotein transmembrane diffusion and clustering events and the iPRD approach is thus a suitable candidate to develop initial models of induced retroviral Env_{3} diffusion.

Here, an iPRD model is developed of Env_{3} diffusion in the context of a truncated planar Gag lattice (Gag_{tr}) confined under lattice-assembly crowding (LAC) conditions. The role of the particular binding stoichiometries and symmetries in the lattice, together with alteration of surface crowding conditions on Env_{3} mobility are considered. The effects of adding reaction-capable PR enzymes are also considered and inform the determinants of induced Env_{3} diffusion coupled to polyprotein degradation.

## 2. Model development

An UCG particle model was developed of truncated Gag polyprotein (Gag_{tr}), including the globular domain of MA and the N-terminal domain of CA (CA_{NTD}) separated by a cleavage site (CS) containing a flexible linker (whose length is based on observed separation of the flanking domains [3,4]), full-length Env_{3} and dimeric PR, using structural dimensions of each of the domains, when known, to inform the radii of corresponding mapped particles (table 1, figure 2*a*).

Effective three-dimensional diffusion constants, , were derived based on a number of considerations including the required dimensionality of diffusion, particle connectivity and assembly conditions (table 1). The diffusion constant for the freely diffusing particle P10 was derived using a modified three-dimensional Stokes–Einstein relation, , where *R*_{P} is the particle radius, *ξ* the dynamic viscosity of water and where the modulation factor *γ*=0.01 was applied to enable a larger simulation time step. It was required for several particles to diffuse two-dimensionally in or at the surface of a membrane (particles P1 and P6–P9). In the iPRD approach, two-dimensional diffusion is achieved by assigning an effective three-dimensional diffusion constant, , to a particle moving under BD, that when constrained by an attractive planar potential (see below), results in two-dimensional motion with an effective diffusion constant, *D**_{2D}. The value of *D**_{3D} for Env_{3} (P6–P9) particles was chosen and tested by preliminary simulations to yield a value for *D**_{2D} compatible with experimental observations [48]. For matrix particle P1, no specific measurements were available. Therefore, parameters for *D**_{3D} were chosen that resulted in values compatible with general experimental ranges for two-dimensional membrane protein diffusion [47]. Particles that were always chemically connected ([P1–P2], [P4–P5,P11], [P6–P9]) were assigned matching diffusion constants. As the CS was always chemically connected until elimination, it was assigned the same value as P1. To emulate the fact that the immature capsid lattice is effectively locked, capsid-related particles (P4, P5 and P11) were assigned diffusion constants, *D**_{3D}, several orders of magnitude lower than that for a single capsid particle in three dimensions, while geometry potentials ensured planarity and spatial confinement (see below). With this choice of diffusion parameters, modulated P10 diffusion still strongly dominates that of other particles involved during the reaction. Thus, it was not necessary to apply modulation factors to other particles.

Space exclusion was applied to all particle types using a harmonic pairwise repulsive (HR) potential (PI1) with force constant *k*_{HR}, centred on the sum of the corresponding pairwise interaction radii (*r*_{I,ij}=*R*_{I,i}+*R*_{I,j}). The form of the HR potential is given by
2.1where *i* and *j* are the two interacting particles, **r**_{i}(*x*_{i},*y*_{i},*z*_{i}) and **r**_{j}(*x*_{j},*y*_{j},*z*_{j}) are their corresponding positions (table 2). Clustering was applied using a weak clustering (WC) potential for inter-P5 (PI2) and inter-P6 interactions (PI3) to model intra-CA hexamer and Env_{3} oligomerization, respectively (figure 2*b*, table 2). The form of the WC potential is given by
2.2where *k*_{WC} is the force constant, *r*_{P,ij} is the sum of the corresponding particle radii and *r*_{C,ij} is the interaction potential cut-off distance.

Rigid oligomerization of MA trimers, MA_{3} and CA hexamers, CA_{6} as well as chemical bonding of disparate domains in Gag_{tr} and Env_{3} were enabled by the assignment of group interactions all consisting of strong pairwise harmonic spring (HS) potentials centred on an equilibrium separation, *r*_{0,ij}, with force constant, *k*_{HS} (figure 2*c* and table 2). The HS interaction potential is given by
2.3

The MA_{3} group (G1), composed of 3 P1, consisted of a single inter-particle interaction between each particle pair in the group. The CA_{6} group (G2), composed of 6 P5, consisted of pairwise interactions between each particle and, clockwise, its: nearest neighbour (GI21), second neighbour (GI22) and furthest neighbour (GI23). In addition, a space-excluding dummy particle, P11, was rigidly assigned to the centre of the CA_{6} hexamer with a HS potential (GI24) between it and each of the P5 particles in the group. The Gag_{tr} and Env_{3} groups (G3 and G4, respectively) consisted of HS potentials between each chemically juxtaposed particle in the group, GI31–GI34 and GI41–GI43, respectively (table 2).

Several particle types were confined by geometry potentials (figure 2*d* and table 2). All potentials were either planar discs centred at coordinates (0,0,*z*_{d}), where *z*_{d} corresponds to the height of the disc above the origin, or cylinders with normal parallel to the *z*-axis and centred at the origin. The centre of initially deposited P1 particles was chosen as *z*=0.

A thick impermeable planar membrane was imposed by two spatially separated repulsive disc (RD) potentials, CI2 and CI5, with force constant, *k*_{RD}, and form:
2.4Calculated placement of particles ensured that no particle initially above the membrane could permeate downwards and vice versa. Planar confinement was imposed on various particles by attractive disc (AD) potentials (CI1, CI3, CI4, CI6 and CI7) with force constant *k*_{AD} and form:
2.5

This combination of potentials ensured that Env_{3} was modelled as a UCG transmembrane protein where neither the ED nor the CT could permeate their respective membrane surfaces while the TMD remained in the membrane. It also enabled Gag_{tr} to be modelled as a monotopic protein representing the tethering by myristoylated MA and explicitly enforced a planar lattice at the level of CA. Boundaries were imposed by repulsive cylindrical (RC) potentials denoted BI1 and BI2 aligned to the *z*-axis, centred at *z*=0, with radius *r*_{B}, height Δ*z*, with force constant *k*_{RC} and with a potential form described by
2.6Different cylindrical boundary radii, *r*_{B}, were used in Study 1 and Study 2, reflecting the different size of the lattice used in each case (table 3). A CA lattice (CA_{lat}) could then be modelled by combining intra-hexamer interactions described by the CA_{6} group with inter-hexamer interactions described by particle interaction potential PI2 and geometry potential CI7, which explicitly enforced a planar lattice geometry on CA_{lat}, confined by cylindrical potential BI1.

Polyprotein decomposition consisted of a bimolecular reaction involving P10 and P3, resulting in P3 elimination (figure 2*e*). As there were no interaction potentials between P2 and P4, elimination enabled disconnect between the P1–P2 fragment and the P4–P5 fragment, corresponding to chemical cleavage of the polyprotein. The experimentally measured rate constant for MA–CA octapeptide cleavage [50] is *k*^{mac}∼*k*_{cat}/*K*_{M}=6.2×10^{3} M^{−1}s^{−1}, corresponding to a value of *k*^{mic}=1.58×10^{2} s^{−1}, in the reaction-limited regime. At physiological concentrations ([MA–CA]=3.8×10^{−3} M, [PR]=0.095×10^{−3} M), a Michaelis–Menten heteropolymer cleavage computation using our earlier formalism [38] yielded 99% completion in approximately 9.6 s, beyond the current computational feasibility for this study. Therefore, a modulated value of *k*^{mic} was assigned, *k*^{mic}=10^{7} s^{−1}, to ensure observable cleavage events within the millisecond time frame of the simulations, yielding *k*^{mac}=3.17×10^{7} M^{−1}s^{−1}, in the diffusion-limited regime.

Initial configurations consisted of deposition of Gag_{tr} chains aligned parallel to the *z*-axis and arranged in a hexameric lattice built around a central hexamer at the origin of the *x*–*y* plane. The MA_{3} symmetry axis lies at the interstice of three different CA_{6}. It is thus not possible to build lattices that simultaneously contain integer numbers of both MA_{3} and CA_{6}. Therefore, two Gag_{tr} lattice types were developed containing: (i) integer CA_{6} with residual MA monomers and (ii) integer MA_{3} with residual CA monomers (figure 2*f*). The *x*–*y* positioning of Env_{3} molecules matched the centre of CA_{6}. Lattice configurations were denoted as [H,R], where H is the lattice ring order (0, single central hexamer; 1, one ring of hexamers around a central hexamer and so on) and R is either 0 for integer numbers of CA_{6}, or 1 depending on if residual Gag_{tr} chains were added to make integer numbers of MA_{3}.

Two studies were performed (Study 1 and Study 2). In Study 1, the effect of imposing MA_{3} and/or CA_{lat} interactions on the diffusion of Env_{3} in the presence of a [1, 0]-Gag_{tr} lattice (systems 1–4) was investigated. In Study 2, the diffusion of multiple Env_{3} was investigated in the absence of Gag_{tr} (system 5), in the presence of a [2, 1]-Gag_{tr} lattice with MA_{3} and CA_{lat} interactions turned on (system 6) and in the additional presence of PR molecules (system 7). Different boundary radii were used in the two studies corresponding to the different lattice sizes (table 3). The boundary radii for systems 1–7 were chosen so that the ratio of the two cylindrical surface areas was maintained at approximately 0.6, corresponding to experimental values, where the inner radius corresponded to LAC conditions and the outer radius to the overall membrane surface. Finally, in Study 2, Env_{3} was investigated in the presence of a [2, 1]-MA_{3} lattice (system 8) under LAC conditions (*r*_{B}(BI2)=33 nm).

The ReaDDY software v. 1.1 was used to implement all iPRD simulations [44]. A time step of 1 ns and temperature of 300 K were used throughout. An ensemble of simulations was performed for each of the systems in Study 1 and Study 2 (table 3). Mean-squared-distances (MSDs) and copy numbers of particles were computed every 1 μs and the time-evolution of the ensemble mean and standard deviations were determined. All simulations were performed on a local compute server; individual ensemble members were simulated on single CPUs, ensembles were computed using a parallelized task-farming protocol.

## 3. Results

### (a) Diffusional locking of Env_{3} by Gag_{tr} lattice

The effect of switching on or off MA_{3} and CA_{lat} interactions on the diffusion of a single Env_{3} initially at the centre of a [1,0] Gag_{tr} lattice was investigated. Switching off MA_{3} interactions signified that G1 interactions were not applied, but MA molecules were still attached to the membrane surface. Switching off CA_{lat} interactions signified that neither CA_{6} (G2) interactions, PI2 interactions, nor CI7 interactions were applied. An ensemble of simulations was performed (figure 3, table 3) for each of four systems corresponding to different permutations of interaction potentials (system 1: [MA_{3} off, CA_{lat} off], system 2: [MA_{3} on, CA_{lat} off], system 3: [MA_{3} off, CA_{lat} on], system 4: [MA_{3} on, CA_{lat} on]).

Switching off both MA_{3} and CA_{lat} interactions (system 1) results in diffusion of Env_{3} with an effective diffusion constant of *D*∼1.5×10^{4} nm^{2} s^{−1} with significant fluctuations (shaded areas represent standard errors). Moreover, even though the CA lattice is confined to a small cylinder (BI1) with radius 18.0 nm, lack of CA_{lat} interactions leads to rapid decay of the hexamer form. Lack of MA_{3} interactions causes MA particles to diffuse as monomers tethered to a confined but disordered CA lattice.

Switching on just MA_{3} interactions (system 2—green) almost completely prohibits Env_{3} diffusion. A [1,0] Gag_{tr} lattice contains only 6 hexamerically arranged MA_{3} molecules (MA_{3} cage) which surround the central Env_{3} molecule. The remaining Gag_{tr} chains have residual MA monomers (not involved in MA_{3} interactions). Even though the CA lattice becomes disordered, confinement of the CA lattice in combination with imposed MA_{3} interactions prohibits the diffusional re-ordering required for particle P9 of Env_{3} to escape the central MA_{3} cage.

Similarly switching on just CA_{lat} interactions (system 3—purple) also almost completely inhibits Env_{3} diffusion. Here, the CA lattice remains hexamerically ordered throughout the simulations with no diffusion as well as being confined. Moreover, MA particles exhibit rapid local fluctuations leading to disordering but being limited by the tethering to the rigid CA lattice. Interestingly, P9 is observed not to escape the central MA cage, even though there are no MA_{3} interactions, suggesting an indirect CA-mediated confinement effect.

Finally, switching on both MA_{3} and CA_{lat} interactions (system 4—black) leads to complete prohibition of Env_{3} diffusion. In this case, induced confinement of MA due to the rigid CA lattice completely limits diffusion of the MA_{3} cage that surrounds Env_{3}. In both system 3 and system 4, the MA layer maintains an approximate although fluctuating hexamer structure.

Therefore, spatial confinement of the CA lattice, its hexameric structure and the trimericity of MA all cooperate to impose rigid diffusional locking of Env_{3} in the membrane. Interaction-absent LAC by Gag_{tr} alone is not sufficient to impede Env_{3} diffusion. Imposing CA_{lat} or MA_{3} interactions is required to severely limit Env_{3} diffusion albeit with local fluctuations, and suggests an indirect stabilizing effect by CA_{lat} that is enhanced by interstitial MA_{3} interactions, while the application of both is completely prohibitive to diffusion.

### (b) Reaction–diffusion unlocking of Env_{3} by PR

Study 2 investigated the effect on Env_{3} diffusion of adding PR enzymes, and thus polyprotein cleavage reactions, to a Gag_{tr} lattice with both MA_{3} and CA_{lat} interactions switched on (figure 4). Simulations were performed for a control system containing 7 Env_{3} and no Gag_{tr} lattice (system 5), a system with 7 Env_{3} and a [2, 1] Gag_{tr} lattice (system 6) and a system with 30 PR additionally added (system 7). Simulations were also performed for a system containing 7 Env_{3} and a [2, 1] MA_{3} lattice (system 8) confined to LAC conditions (*r*_{BI2}=33 nm) to compare with the effect of reduced crowding.

In system 5, Env_{3} diffuses with *D*∼1.25×10^{4} nm^{2} s^{−1} (green) similar to the value in figure 3 when both MA_{3} and CA_{lat} interactions are switched off. Furthermore, clustering of at least 2 or 3 Env_{3} is exhibited in most simulations within 5 ms (electronic supplementary material, movie M1). Diffusion of Env_{3} in system 6 is completely prohibited (black). Although each Env_{3} maintains rapid local fluctuations, no Env_{3} molecules are able to escape the intra-hexamer holes in which they are originally placed. Furthermore, MA particles are indirectly locked within a radius of 33 nm, imposed by the explicit confinement of the CA lattice, in all simulations (electronic supplementary material, movie M2).

Addition of PR (system 7) results in a switch from an initially diffusion prohibited system to one which diffuses with *D*∼8.33×10^{3} nm^{2} s^{−1}, only approximately 1.5-fold slower than with Env_{3} molecules alone. The diffusion switch occurs within the first 1 ms of simulations (figure 4*a*) and is concomitant with the decay of the Gag_{tr} upon CS elimination by PR (figure 4*b*). It is also concomitant with the diffusional spread of MA_{3} molecules decoupled from the CA lattice to fill the larger bounded surface area, with radius 42 nm. This results in decreased crowding as followed by the change in MA inner-surface density from *ρ*=0.04 nm^{−2} to *ρ*=0.03 nm^{−2} (electronic supplementary material, movie M3). Cleavage of Gag_{tr} at the MA–CA junction thus enables decoupling of the MA layer from the CA lattice, switching on Env_{3} diffusion. Env_{3} diffuses with *D*∼5×10^{3} nm^{2} s^{−1} (blue) in a LAC-confined MA_{3} layer (system 8), corresponding to *ρ*=0.04 nm^{−2}. Therefore, the effect of reduced crowding as a result of lattice decoupling contributes a 1.7-fold speed-up in diffusion compared with decoupled LAC confinement alone.

## 4. Discussion and conclusion

Using HIV-1 as a model system, the basic reaction–diffusion mechanisms at work during the induction of retroviral infectivity have been investigated here by the development of a mesoscopic spatio-temporal computer model of lateral membrane diffusion by Env trimers (Env_{3}) in the presence of a monotopic polyprotein lattice (Gag_{tr}) composed of matrix–capsid (MA–CA) protein domains that is chemically decomposed by protease (PR) enzymes.

The results reported here suggest that, in the absence of both trimeric MA interactions and CA lattice interactions, maintaining LAC conditions by Gag_{tr} is not alone sufficient to significantly impair Env_{3} diffusion. Diffusional locking of Env_{3} occurs substantially in LAC conditions with Gag_{tr} if either MA can form trimers (MA_{3}) or CA can form a hexameric lattice (CA_{lat}).

LAC conditions in combination with MA_{3} and CA_{lat} interactions representing those that exist in a retroviral Gag lattice, result in complete locking of Env_{3} within the MA layer. This suggests MA and CA multimericities cooperate to rigidly lock Env_{3} in the membrane. The CA lattice induces an N-ward effect on MA intra-hexamer rigidity and thus indirectly on Env_{3} mobility, while MA trimerization induces a C-ward effect on CA, enhancing inter-hexamer stability. This could be tested by Env_{3} mobility experiments under mutationally induced variations of MA and CA multimer binding affinities.

Although MA forms hexamers of trimers within *in vitro* lattice arrays, neither intra-hexamer MA interactions nor long range order in the MA layer has been identified *in vivo* [7]. It is suggested here that confinement in initial hexamer of trimer form may exist *in vivo* but is mainly consequential on the underlying immobile CA lattice. Adding PR results in cleavage of the inter-MA–CA linker within Gag_{tr}, decoupling diffusion of MA from CA. The results reported here show that, in the absence of long range matrix order, this decoupling enables MA_{3} diffusion in turn inducing Env_{3} diffusion. The effect of decoupling alone—calculated by Env_{3} diffusion under maintained LAC conditions of MA_{3}—is sufficient to recuperate diffusivity to 40% of MA_{3}-absent Env_{3} diffusion. The combination of decoupling and consequential decrowding recuperates Env_{3} diffusivity to 66% of the MA_{3}-absent speed. Therefore, Gag_{tr}-decoupled MA_{3} molecules at physiological concentrations only partially interfere with Env_{3} diffusion.

In order to make the simulations computationally feasible, accelerated reaction rate constants were used in this study. It is unlikely that altering reaction rate would qualitatively change the striking differential effect observed between a PR-containing and PR-absent system—that is Env_{3} does not diffuse in the PR-absent system. Thus the PR reaction is essential to induce a diffusional switch in Env_{3}.

However, the quantitative time scale would be dramatically affected and this may change the nature of the diffusional induction. In the accelerated simulations, the reaction completes (within 1 ms) long before Env_{3} can form even incomplete clusters (within 5 ms). Thus, the diffusional switch is sharply reaction-mediated. Using measured *in vitro* cleavage rate constants, Gag_{tr} processing via Michaelis–Menten degradation under *in virio* concentrations would take approximately 10 s. Processing the whole of Gag would take considerably longer, owing to the differential and often slower rates for various cleavage junctions [38]. Total processing of the MA–CA junction during autocatalysed Gag and GagPol degradation under similar conditions has been computed to take as long as 20 min [39]. While crowding and architectural conditions in the virion may accelerate this, it is unlikely to be faster than the seconds time scale.

Unimpeded, Env_{3} molecules that are diffusionally liberated (*D*_{2D}∼1.25×10^{4} nm^{2} s^{−1}) during the reaction process could potentially explore the virion surface within approximately 5 s, and may therefore be able to encounter still-constrained Env_{3} molecules. Diffusion of liberated Env_{3}, however, would be modulated by limited surface area access for some time into the reaction, owing to the still existent Gag lattice. This would in turn limit access to still-constrained Env_{3}. Nonetheless, this suggests that Env_{3} oligomerization completes rapidly during the PR reaction once both sufficient Env_{3} are liberated and a significant degree of MA–CA lattice decoupling has taken place. This lends to a more broadly reaction-mediated process whose controlling parameters could be explored further.

Spreading of MA_{3} molecules across the entire membrane area results in reduced crowding, which slightly enhances Env_{3} diffusion (1.7-fold calculated here) compared to that in LAC conditions of MA_{3}. This initially suggests the role of partial Gag coverage observed in cryo-EM experiments of budded virions [8] could be to provide diffusional control of Env_{3}. However, at physiological reaction time scales of minutes a slight speed-up in Env_{3} diffusion and thus surface exploration on the time scale of seconds would not impact the overall time to achieve clustering. This rather suggests that partial Gag coverage plays a limited role in diffusional Env_{3} control. Nonetheless, experiments in which Env_{3} diffusivity is measured at various lattice coverage ratios could test whether either Env_{3} diffusion or the time to reach infectious Env_{3}-oligomerized virions is significantly modulated.

Furthermore, a number of assumptions have been made on the level of coarse-graining, the planar geometry, the radii of particles, the strength of interaction potentials and the effective particle diffusion constants assigned. While these assumptions have been informed by available data, future work could parametrize these values against emergent diffusion experiments.

Despite its simplicity, the Env_{3} diffusional switch model purported here explains the observed coupling between PR-dependent Gag degradation and induction of fusion activity [28], from an underlying physical perspective. The model is consistent both with observations that truncation of the Env_{3} cytoplasmic tail (CT) can induce fusion activity [24] and that Env_{3} molecules cluster on the surface of mature but not immature virions [29]. Effective CT size may then be of considerable interest for regulating infectivity. The CT structure is currently unknown. Interestingly, its sequential size varies significantly across different retrovirus species [25]. If sequence length is indicative of spatial extent, then it would imply that retroviruses have evolved CTs to regulate Env_{3} diffusion induction and thus the time scale to the onset of virion infectivity.

From a synthetic retrovirology perspective, understanding the diffusional and aggregational determinants of Env_{3} in terms of a physical modelling framework may assist in the design of novel pseudotyped retroviruses and virus-like particles (VLPs) for use as gene delivery vectors [51] and vaccines [52]. Moreover, accelerated reaction conditions could be engineered through mutatgenesis (*k*_{cat} enhancement to 10^{4}–10^{5} s^{−1} would still be well within enzymatic ranges) to a regime where diffusional induction was sharply reaction-mediated. Then slight variations in diffusional properties, surface coverage and reaction rates could sensitively control infectivity switching. In principle, logic gates could be constructed from VLPs and suitable input-signal molecules. For example, an Env_{3}-separated Gag_{tr}-locked retroviral VLP (output signal: 0), with viral PR and a PR inhibitor serving as molecular input signals, is conceptually equivalent to an INH logic gate. Only input of PR and not its inhibitor (input signal: 1 0) would result in reaction-mediated diffusional unlocking and clustering of Env_{3} in the VLP (output signal: 1).

The reaction–diffusion model developed here is extensible. Development of multiparticle capsid interaction potentials is forthcoming, the aim being to model capsid self-assembly with sufficient accuracy within an iPRD approach. This would enable pathways of capsid self-assembly during the proteolytic reaction to be followed and may provide additional insight into the assembly mechanism.

The full process of retroviral maturation remains a formidable challenge to simulate. The systems studied here are deliberately simplified truncated systems, but still reach of the order of 1000 particles and achieve ensembles of milliseconds on local CPU clusters. At this level of UCG, an entire virion may require of the order of 10^{4}–10^{6} particles, while maturation time scales are of the order of minutes. Recent advances in GPU implementations of iPRD simulations [53] provide substantial acceleration and, in conjunction with larger compute resources, mesoscopic Markov state models [54] and structural packing methods [55], may make modelling such systems approachable in the near future. However, it is likely that multiparticle potentials as well as hierarchical multiscale modelling approaches that combine iPRD simulations with coarse-grained MD and rigid-body diffusional association of proteins [42] may be required to adequately model phenomena such as the complete process of reaction–diffusion-coupled polyprotein degradation, capsid self-assembly and ribonucleprotein packaging that all take place during virion maturation.

## Competing interests

The author declares that he has no competing interests.

## Funding

The author acknowledges support from amfAR Mathilde Krim Fellowship in Basic Biomedical Research number 108680.

## Acknowledgements

The author is grateful to Dr Natalia Gabrielli for critically reading the manuscript, Jan Huertas for technical assistance and to Prof. Andreas Meyerhans, Prof. Gilles Mirambeau, Dr Sébastien Lyonnais, Prof. Frank Noé, Dr Johannes Schöneberg, Dr Alexander Ullrich and Johann Biedermann for insightful discussions and suggestions. The author also thanks the referees for their valuable suggestions during the review process.

## Footnotes

One contribution of 17 to a theme issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

- Accepted August 22, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.