To Build Viruses: A Markov Chain Monte Carlo 
Algorithm for Simulating Viral Assembly Kinetics

Brunk, N(SuperDagger[.]) , Butske, W((.)^‡) , Zlotnick, A((.)^(SuperStar[‡]))

†Department of Physics, Indiana University
‡Department of Mathematics, Rose-Hulman Institute of Technology
(SuperStar[‡])Department of Molecular & Cellular Biochemistry, Indiana University

Join the Conversation #WolframTechConf

Physical Virology:  Capsid Self-Assembly Reaction Kinetics

What is a capsid?

What role does it play in the virus life cycle?

Capsid:  a protein shell made of many copies of few (or even singular) types of subunits forming closed polymers with icosahedral or helical symmetry.

Role:  Shell protects fragile viral cargo (genome, transcription enzymes, etc.) during intracellular transfer.

Free energy minimization -> spontaneous self-assembly.

	Example:  Hepatitis B Virus Capsids (120, 90 subunit isomers).	[actual  T = 4 has rotated hexamers, see following slide]

Leftmost image courtesy:  cgl.ucsf.edu		All other images, unless otherwise stated, created per this work.

Reaction Kinetics Mathematics

A Toy Example:  Tetrahedron Assembly

Simplest case: self-assembly of a tetrahedron.  Sequentially adding a single subunit at a time:

Mass conservation -> (Subscript[z, 1])' == - (2(Subscript[z, 2])’ + 3 (Subscript[z, 3])’ + 4 (Subscript[z, 4])’) , thus forming a 4 state variable ODE system.

Definitions:  Forward Reaction Rate, Statistical Terms

Forward reaction rate ((Subscript[k, f])) is “microscopic;” statistical terms account for the macroscopic system behavior in vitro.

Statistical terms {δ, (Subscript[s, i]), (Subscript[σ, i])} account for the combinatorics of subunits, scaling all rates linearly:

(δ) is the number of binding sites on fundamental subunit ((Subscript[z, 1])),

((Subscript[s, i])) is the number of sites upon which binding forms intermediate ((Subscript[z, i])) from ((Subscript[z, i-1])),

((Subscript[σ, i])) is the number of bound subunits in ((Subscript[z, i+1])) whose dissociation forms ((Subscript[z, i])).

Definitions:  Backwards Reaction Rate, Energetics

Backward reaction rates ((Subscript[k, (Subscript[b, i])])) := ((((Subscript[k, (Subscript[b, 0])])))/((Subscript[K, (Subscript[A, con])])^(Subscript[χ, i]))) , where single-contact association constant (Subscript[K, (Subscript[A, con])]) = (e^(-(((((Subscript[ΔG, con])/((Subscript[k, B]) T))))^β))) , a stretched Boltzmann Factor / Arrhenius relation.

((Subscript[ΔG, con])) is the energy of a single contact,

((Subscript[k, b]) T) is the standing thermal free energy,

(β) represents a superimposition of slightly variable ((Subscript[ΔG, (Subscript[con, i])])’s),

& ((Subscript[χ, i])) is the number of contacts formed upon creation of (Subscript[z, i]).

Tetrahedron combinatorics are trivial.  However, that of larger structures are not:

Dodecahedral Assembly:  A Simplest Case 12 Equation Model

Tetrahedron above:  a toy used to explain the mathematics in a comprehensible format.  

Viruses form more complex polyhedra, often consisting of hundreds of subunits.

	Bovine Enterovirus forms a dodecahedron of ~ 16nm diameter.

Most-stable subset thereof forms a 12 state variable (equation) model displayed below:

GraphicsRow[autoRotate/@intPolyView/@First/@First/@Uncompress/@bank,ImageSize->Full]

Graph Theory:  Computational Tractability of Assembly

Graph Theory Introduction

What is a graph?

A set of vertices {(Subscript[v, i])} joined by a set edges {(Subscript[e, i])}.

	Examples:	(note Graph3D[~] functionality chronology, Mathematica Stack Exchange)

GraphicsRow[{
Graph[{1\[UndirectedEdge]2,2\[UndirectedEdge]3,3\[UndirectedEdge]1},GraphHighlight -> {1\[UndirectedEdge]3}, VertexLabels->"Name",EdgeStyle->Thick],
Graph[{1\[UndirectedEdge]2,1\[UndirectedEdge]3,1\[UndirectedEdge]4,2\[UndirectedEdge]3,2\[UndirectedEdge]4,3\[UndirectedEdge]4},GraphHighlight -> {1\[UndirectedEdge]3}, VertexLabels->"Name",EdgeStyle->Thick,
VertexCoordinates->PolyhedronData["Tetrahedron","VertexCoordinates"]],
Magnify[VertexDelete[Graph3D[{1\[UndirectedEdge]2,1\[UndirectedEdge]3,1\[UndirectedEdge]4,2\[UndirectedEdge]3,2\[UndirectedEdge]4,3\[UndirectedEdge]4},GraphHighlight -> {1\[UndirectedEdge]3}, VertexLabels->"Name",
VertexCoordinates->PolyhedronData["Tetrahedron","VertexCoordinates"]],{}],1.4225]},
ImageSize->Full]

Capsid Assembly & Graph Theory:  A Connection

Note mathematical property of duality:

Graph theory renders assembly combinatorics computationally tractable.

Capsid subunits:  vertex set {(Subscript[v, i])} ,

Subunit contacts:   edge set {(Subscript[e, i])}  :

EdgeCount[intG] = energy within intermediate.

Rendering contact formation (& relevant statistical parameters) computatioanlly tractable.

	-> energy minimization pathways during assembly.
	
To keep track of either of these (HBV variants), you need a computer:

Enumeration:  Defining State Space

Model State Space (ie. an Intermediate Bank)

With models and in vitro, kinetic trapping occurs in which too high a value for ((Subscript[ΔG, c])) causes aggressive enough assembly to deplete the required subunit pool prior to appreciable capsid assembly, trapping otherwise transient intermediates.  

Prediction of intermediate activity is a primary goal of this work.  Adequate predictions require more thorough physical treatment.  

Excessive assumption:  that only the most stable intermediate of each size contributes to the system.

More realistic system : (73) state variables.

An “assembly funnel:”

State Space:  Alternative Views

There are multiple ways in which one can view the state space.  

An energy diagram representation displaying the relative probabilities of generation (via the assembly algorithm) is displayed below:

(* Plots the intermediates generated:  *)
With[{ξ=N@Table[(# /. {a_List, count_} :> {zi,((-(10^count))/10),-(((EdgeCount[intG[a]])/(Length[a]))/(12/5))}) &/@Uncompress@bank[[zi]],{zi,1,12,1}]},
Show[ListPointPlot3D[ξ,PlotTheme:>{"Marketing","Detailed"},PlotRange:>All,PlotStyle:>PointSize[.0075],AxesLabel:>Evaluate[Style[#,16]&/@{"(Subscript[z, n])","- Log(ϕ)","((Subscript[ΔG, c])/n)"}],ViewPoint->{(Pi/1),((-Pi)/1),(Pi/6)},ImageSize:>1000,PlotLabel->"Dodecahedron Assembly Intermediates (73)"]
(*Graphics3D[{Black,Line[Flatten[ξ,1]]}]*)]]//Rasterize

The Assembly Pathways Graph:  A More Utilitarian View

A graph with: each intermediate constituing a vertex and the possibility of transition between the two (via single subunit addition) constituing an edge.

vertices:  isomerizations possible (above)

edges:  possibility of transition between them by single subunit addition / removal.

The dodecahedron’s pathway graph (omitting central, backwards edges):

Effect on the Math

Magnify[GraphicsRow[{waysUpTableDisp[5],waysDownTableDisp[4]},ImageSize->Full],2]

Magnify[showRHSps[6,5,100],2]/.transform

The Enumeration Algorithm

Functionality, Explanation

A Markov Chain Monte Carlo assembly algorithm allows for enumeration of the above “intermediate bank” ie. the state space of the model.  

Each memoryless, stochastic Monte Carlo step adds a single subunit:

Manipulate[intPolyView[u1[[1;;ν]],{},{ImageSize:>Large}],{ν,1,12,1},Button["New Trajectory",u1=assemble[12]],SaveDefinitions->True]

Isomerization Sieve Rules

(2) Isomerization Rules:

Graph isomorphism:  a graph mapping query pertaining to combinatorics without symmetry.

Yields 57 / 73 intermediates, all energetically distinguishable paths.

Used elsewhere in literature (Hagan, etc.)

Custom sameOrbitQ[(Subscript[τ, 1]),(Subscript[τ, 2])] query.

Yields remaining 26 isomers.

Uses group theory:  are any ((Subscript[τ, 1])) & ((Subscript[τ, 2])) of the same orbit under permutation?

If so, the candidates are the same isomer.

Exhaustive Dodecahedral Model (73 state variables)

Future Directions

State Decomposition Allows Scaling

Virus capsids may consist of hundreds or more subunits; isomer count for a 30-mer is approximately 2.6 million.

Thus, scale-up strategy must involve state decomposition ie. complexity reduction.

In a stability-sorted “assembly funnel,” ((system change)/(state variable added)) yields diminishing return with inclusion depth:

May safely neglect highly unstable majority of intermediates.

Stochasticity of each Monte Carlo step allows preferential selection of subunits that yield more stable isomers.

Thus, conjoining most stable Markov Chains, will provide sparse models approximating larger virus assembly.

	(This will ensure connectivity as well.)

Analytical Characterization & Predictions

Capsid Kinetics & Curve Fitting

In vitro activity may be observed by a variety of experimental techniques, both well established and novel:

	Standard techniques for the observation of capsid kinetics include:

Electrophoresis Gel Densitometry

Size Exclusion Chromatography

Dynamic Light Scattering

These will provide a means of observing and curve-fitting kinetics data for only the capsid activity.

	Scalability allows characterization of different viruses’ kinetics.  For example, a recent virus’ kinetics fit to the corresponding model (species, model size undisclosed):

Demonstrable that this is within the expected max precision of curve fitting partial populations.

Seeking combination of parameters that yields the least variability.

Seeking Correlation w/ Intermediate Activity

Novel techniques for observation of intermediate activity are in development at IU, including:

Charge Detection Mass Spectrometry

Nanofluidics

These will provide a means of seeking intermediate activity that correlates with the mathematical systems.

Intermediate activity observed thus far by multiple labs has correlated well with preliminary models, specifically:

very early time, disproportionate occurrence of less stable intermediates due to recurrent pathway combinatorics,

appreciable buildup of intermediates late in the assembly pathway, ie. near (Subscript[z, capsid]).

Contact information:  Nicholas Brunk, nbrunk@indiana.edu