The output of the simulation can be visualised using various graphs, which are updated in real-time. These plots are displayed in the Plots panel and in the Sitewise panel. The display of each plot can be configured with its respective settings button. The raw data (.tsv) or a high resolution image (.png) can be downloaded . The different plots available are detailed below.

If uncertain about the value of a parameter, it can be given a background / prior distribution to capture its uncertainty. The probability distributions supported are uniform, normal, lognormal, and exponential. The parameter is resampled at the beginning of each simulation. All parameters are sampled independently of each other. This feature is useful for exploring the effect that parameters have on the system, or for fitting to a dataset.

The Gibbs energy $G$ describes the work that can be performed by a thermodynamic system. The standard Gibbs energy of a reaction $\Delta G$ describes the change in Gibbs energy that accompanies a reaction of 1 M of reactants, under standard conditions (at 25°C and 100 kPa). $\Delta G_s$ describes how energetically favourable state $s$ is. A system spends more time in states which have low Gibbs energies. Suppose we have a single molecule which can exist in a number of chemical states. When the system has reached equilibrium, the probability of the molecule being in state $s$ at any given time is estimated with a Boltzmann distribution:

$$ p(s) \propto e^{-\frac{\Delta G_s}{k_B T}}, $$

where $\Delta G_s$ is the Gibbs energy of state $s$, $k_B$ is Boltzmann's constant ($1.3806 \times 10^{-23}$J.K$^{-1}$), and $T$ is the temperature in Kelvins. Energy is expressed in units of $k_B T$, by taking the experimental energy (in units of J) and dividing by $k_B T$.

For example, suppose a single molecule is in equilibrium between two states, $A$ and $B$. This is expressed as $A \rightleftharpoons B$. Suppose that the states have Gibbs energies of $\Delta G_A = -3\;k_B T$ and $\Delta G_B = -6\;k_B T$, respectively. We would therefore expect the molecule to be in state $A$ 4.7% of the time:

$$p(A) = \frac{e^{-\frac{\Delta G_A}{k_B T}}}{e^{-\frac{\Delta G_A}{k_B T}} + e^{-\frac{\Delta G_B}{k_B T}}} = \frac{e^{-\frac{-3\;k_B T}{k_B T}}}{e^{-\frac{-3\;k_B T}{k_B T}} + e^{-\frac{-6\;k_B T}{k_B T}}} = \frac{e^{3}}{e^{3} + e^{6}} = 0.047.$$

The Boltzmann distribution assumes that the system has been fluctuating for long enough to achieve equilibrium. However, since the competing processes of transcription (see Reactions) occur on similar timescales, equilibrium approximations do not always make good assumptions in this context.

Chemical kinetics is the study of rates of chemical processes, and how long reactions take to occur. A chemical reaction can be described with a rate constant $k$ (units s$^{-1}$). $k$ describes the speed of the reaction - a higher rate means that more reactions are expected to occur in a given time. Energy is required for a system to transit from one state to another. The standard Gibbs energy of activation for a reaction is denoted by $\Delta G^\dagger$. The larger this barrier, the more time needed for random thermal collisions to provide enough energy to allow the system to make the transition. Suppose a system is currently in state $s$, with Gibbs energy $\Delta G_s$. The Arrhenius equation can estimate the rate constant of transiting from state $s$ to state $s^\prime$:

$$ k_{s \rightarrow s^\prime} = A e^{-\frac{\Delta G^\dagger - \Delta G_s}{k_B T}}, $$

for some pre-exponential constant (units s$^{-1}$), where $\Delta G^\dagger$ is the absolute height of the Gibbs energy barrier between states $s$ and $s^\prime$. Once the rate constants of all possible reactions are known, we can start the single-molecule simulation.

Transcription is simulated using the Gillespie algorithm. A critical assumption of this method is that these reactions are Poisson processes; the amount of time until a reaction occurs is independent of all other reactions and independent of the amount of time elapsed. For example, suppose we were to continue flipping a coin until it lands on tails. The probability of getting tails on the next flip is not affected by how many times we have already flipped, and is not affected by how many times someone else has flipped their coin. Similarly, whether a chemical reaction has been waiting 1 second or 1 hour for the reaction to occur, the probability of it happening right now is always the same. The amount of time until reaction $s \rightarrow s^\prime$ occurs follows the exponential distribution with rate $k$.

If there are multiple chemical reactions competing with each other, then the reaction with the shortest reaction time will be selected. For example, suppose we have have a test-tube containing a single molecule of compound $A$, which can spontaneously undergo two reactions: $A \rightarrow B$ or $A \rightarrow C$, with rate constants 0.001s^-1 and 0.002s^-1 respectively. The test-tube will flash a different colour of light for each reaction.

Adding the next base onto the $3^\prime$ end of the nascent strand requires one translocation step (see Kinetics of translocation) and at least three chemical steps:





Any of these steps could be reversible. As there is not sufficent information to model all three steps, nucleotide incorporation has been compressed into a two step process.




(units s$^{-1}$μM$^{-1}$) is the second order rate constant of binding NTP. The rate of binding an NTP is equal to $\times$ . Once the NTP has been bound, it may be released back into solution. The rate of release is calculated from the dissociation constant (units μM). Rate of NTP release is equal to $\times$ .

When the NTP concentration is high, or if is large, or if NTP binding is approximated as an equilibrium process ($ = 1$), then NTP binding becomes so rapid that equilibrium between these two states is achieved. Under these conditions, the kinetics of NTP binding are no longer relevant. is simply the ratio between the rates of binding and release, and can be used to calculate the proportion of time that NTP is bound. For example, suppose that = 35μM and = 1000μM, then the proportion of time that NTP is bound (while in this two-state system) is around 97%:

$$ p(\text{bound}) = \frac{ \frac{\text{[NTP]}}{K_D} }{\frac{\text{[NTP]}}{K_D} + 1} = \frac{ \frac{1000\text{ μM}}{35 \text{ μM}} }{\frac{1000\text{ μM}}{35 \text{ μM}} + 1} = 0.966. $$

Approximating NTP binding/release as an equilibrium process (using ) is equivalent to setting to a high value, however is much faster to simulate.

Let $\Delta G_{S}$ be the Gibbs energy of state $S$. In the context of nucleic acids, this can be approximated using nearest neighbour models. The basepairing energetic contribution of a transcription state is computed by summing the Gibbs energy of the gene to that of the hybrid:

$$ \Delta G_{total}^{(bp)} = \Delta G_{gene}^{(bp)} + \Delta G_{hybrid}^{(bp)} $$

Let $\Delta G_{S(l, t)}$ be the Gibbs energy of state $S(L, t)$. This term is calculated using:

$$ \Delta G_{S(l, t)} = \Delta G_{S(l, t)}^{(bp)} + \begin{cases} \Delta G_{\tau 1} &\text{ if } t = 1 \text { (posttranslocated) } \\ 0 &\text { if } t \neq 1 \end{cases} $$

where $\Delta G_{S(l, t)}^{(bp)}$ is computed from themodynamic parameters of nucleotide basepairing and is the Gibbs energy added onto the posttranslocated position.

If the rate of translocation is very rapid, or if translocation between the pre and posttranslocated states is approximated as an equilibrium process ($ = 1$), then these energetic terms describe the proportion of time the system spends in the two states. This is described by translocation equilibrium constant $K_\tau$.

$$ K_\tau = \frac{k_{S(l,1) \rightarrow S(l,0)}}{k_{S(l,0) \rightarrow S(l,1)}} = \frac{k_{bck}(S(l,1))}{k_{fwd}(S(l,0))} $$

However, in the general case, translocation is not assumed to reach equilibrium ($ = 0$). Basepairs may need to be broken in order for the polymerase to translocate forwards (right) or backwards (left) from its current position. For example, consider the state below.




Let $S(l,t)$ and $S(l,t+1)$ be two neighbouring states and let $T(l,t)$ be the translocation transition state between the two. Using transition state theory, the Gibbs energy of this transition state is calculated as:

$$ \Delta G_{T(l, t)} = \Delta G_{T(l, t)}^{(bp)} + \Delta G_\tau^\dagger + \begin{cases} \Delta G_{\tau -}^\dagger &\text{ if } t \leq -1 \text{ (backtracking) }\\ 0 & \text { if } t = 0 \text{ (elongation) } \\ \Delta G_{\tau +}^\dagger &\text { if } t \geq 1 \text{ (hypertranslocation) } \end{cases} $$

where $\Delta G_{T(l, t)}^{(bp)}$ is the Gibbs energy of basepairing in the transition state, is the Gibbs energy barrier of translocation, is the Gibbs energy barrier of backtracking, and is the Gibbs energy barrier of hypertranslocation. Estimating $\Delta G_{T(l, t)}^{(bp)}$ can be achieved through several methods (discussed below). After the translocation energy barrier has been estimated, the rates of translocating between $S(l,t)$ and $S(l,t+1)$ can be calculated using the Arrhenius equation:

$$ k_{S(l,t) \rightarrow S(l,t+1)} = k_{fwd}(S(l,t)) = A e^{-(\Delta G_{T(l, t)} - \Delta G_{S(l, t)})} \\ k_{S(l,t+1) \rightarrow S(l,t)} = k_{bck}(S(l,t+1)) = A e^{-(\Delta G_{T(l, t)} - \Delta G_{S(l, t+1)})} $$

using pre-exponential constant $A = 10^6$ s^-1.

SimPol has several methods for estimating $\Delta G_{T(l, t)}^{(bp)}$ (models are denoted by ):



In optical trapping experiments, beads are tethered to single molecules of RNA polymerase and optical tweezers are used to exert a force (units pN) upon the elongation complex. Positive forces $= F > 0$ act in the same direction as transcription (an assisting load) while negative forces $F < 0$ act in the opposite direction (a hindering load). Increasing the assisting force serves to decrease the Gibbs energy barrier of forward translocation and increase the barrier of backwards translocation. This process is modelled by adding an additional term to the Arrhenius equation:

$$ k_{fwd}(S(l,t)) = A e^{-(\Delta G_{T(l, t)} - \Delta G_{S(l, t)} - F(\delta_1))} \\ k_{bck}(S(l,t)) = A e^{-(\Delta G_{T(l, t-1)} - \Delta G_{S(l, t)} - F(\delta - \delta_1))} $$

where $\delta$ is the distance between consecutive basepairs ($\delta = 3.4$ Å in dsDNA), and 0 < < $\delta$ is the position of the transition state between the two neighbouring states.

The translocation equilibrium approximation asserts that is so small that translocation between the pretranslocated ($t=0$) and posttranslocated ($t=1$) states occurs on a timescale significantly faster than other reactions eg. NTP binding and incorporation.

Let $D$ be a dataset and let $\Theta = ${$\theta_1, \theta_2, \dotso, \theta_b$} be a set of parameters for a model. The goal of statistical inference is to find values for each parameter in $\Theta$ such that the model adequately describes the data $D$.

In Bayesian statistics, parameters are estimated using the posterior distribution. The posterior probability of $\Theta = \Theta_k$ can be obtained using: $$P(\Theta_k|D) = \frac{P(D|\Theta_k)P(\Theta_k)}{P(D)}$$ where $P(D|\Theta_k)$ is the likelihood - the probability of observing data $D$ under the model assuming that the parameters in $\Theta_k$ are set to their true values - and $P(\Theta_k)$ is the prior probability of $\Theta_k$ being the true values before observing new data $D$. $P(D)$ is a normalisation constant and is often intractable to evaluate directly. Subsequently, sampling methods are typically used to estimate the posterior distribution $P(\Theta|D)$. In a Bayesian framework, parameters are estimated as probability distributions - ie. the posterior distribution - instead of as point estimates.

Unfortunately, in some circumstances calculating the likelihood $P(D|\Theta)$ is not tractable. Instead, simulation can be used, and the posterior distribution can approximated. This framework is known as approximate Bayesian computation (ABC). Let $\rho(D,\hat{D}_k)$ be a function which evaluates the distance between the true data $D$ and the data $\hat{D}_k$ simulated under the model with parameters $\Theta_k$. Then, $\Theta = \Theta_k$ is accepted into the posterior distribution if $\rho(D,\hat{D}_k) \leq \epsilon$ for some threshold $\epsilon \geq 0$. The posterior distribution approximated under $\epsilon$ is denoted by $\hat{P}_\epsilon(\Theta_k|D)$.

If $\epsilon$ is sufficiently small, then $\hat{P}_\epsilon(\Theta_k|D)$ will be an accurate approximation of the true posterior distribution. However, if $\epsilon$ is too small, then adequately sampling from $\hat{P}_\epsilon(\Theta_k|D)$ may not be attainable in any reasonable timeframe.

In this section, first the various types of Experimental Data $D$ supported by SimPol are described. Then, two ABC algorithms for estimating $\hat{P}_\epsilon(\Theta_k|D)$ are explained: R-ABC and MCMC-ABC.

SimPol supports three types of data, detailed below. The total experimental data is comprised of a series of $c$ experiments: $D = (d_1, d_2, \dotso, d_c)$. In the image above, there are $c=2$ experiments. Each datatype is associated with a distance function $\rho$. The values of $\rho$ (rho) are summed across each of the $c$ experiments to give the total distance.

The rejection-ABC (R-ABC) algorithm is simple and easy to parallelise. Each sample is generated independently. The following user-specified settings are involved:

• N iterations, $N$: the total number of samples from posterior distribution $\hat{P}_\epsilon(\Theta|D)$ to generate.
• N trials per data point, $n$: the number of simulations per data point. The average measured value (mean velocity or median pause time) will be taken.
• $\epsilon$: the threshold that $\rho(D,\hat{D}_k)$ must meet in order for $\Theta_k$ to be accepted into $\hat{P}_\epsilon(\Theta|D)$. This setting is linked with and can be used instead of $P(\rho \leq \epsilon)$.
• $P(\rho \leq \epsilon)$: the probability that a value sampled from $P(\Theta)$ will be accepted. This setting is linked with and can be used instead of $\epsilon$.

To test this algorithm on motivating example 1, follow the link below and then find and press . The posterior distribution can be summarised into a single state using .

R-ABC motivating example 1 session: Open | View XML

The Markov chain Monte Carlo ABC (MCMC-ABC) algorithm produces a random walk through parameter space. Samples are no longer independent and therefore autocorrelation between successive states should be accounted for. In contrast to R-ABC, $\epsilon$ is not fixed. To enable $\rho$ to converge, $\epsilon$ starts out high at $\epsilon_0$ and then decreases to $\epsilon_{min}$ following an exponential cooling scheme: $\epsilon_i = \max(\epsilon_{min}, \gamma \epsilon_{i-1})$, where $0$ < $\gamma$ < $1$. $\epsilon_i$ refers to the value of the threshold at iteration $i$ in the MCMC chain. If $\gamma$ is too small, the system may not converge, whereas if $\gamma$ is too large then it may take too long to converge. The following user-specified settings are involved:

• N iterations, $N$: the length of the chain of posterior distribution samples $\hat{P}_\epsilon(\Theta|D)$.
• N trials per data point, $n$: the number of simulations per data point. The average measured value (mean velocity or median pause time) will be taken.
• $\epsilon_{min}$: the final value of $\epsilon$ after cooling.
• $\epsilon_{0}$: the initial value of $\epsilon$.
• $\gamma$: the rate of $\epsilon$ cooling.
• Burn-in (%): the percentage of samples to discard before reaching the posterior distribution. Set to -1 to allow automatic detection.
• Log every (states), $L$: the frequency of storing posterior samples. Storing every sample is wasteful because of autocorrelation between successive samples.

To test this algorithm on motivating example 1, follow the link below and then find and press . The posterior distribution can be summarised into a single state using . The chain is not guaranteed to converge and may tay a long time even if it does.

MCMC-ABC motivating example 1 session: Open | View XML

Recommended only for advanced users: the model itself can be estimated in Approximate Bayesian computation as well as the various Parameters. Kinetic models, which may differ in their model settings (for example equilibrium assumptions, whether backtracking is enabled) or by constraining parameters at fixed values, can be specified in the Model comparison panel. Prior probabilities may be assigned to each model by giving each model a weight. The probability that the model will be activated at the beginning of each simulation, or the prior probability of that model during approximate Bayesian computation, is proportional to that weight.

Upon pressing a new kinetic model will be created. The settings of that model will be those that are currently specified. Any parameters which do not have a prior distribution will be held fixed at their current value, and those which do have a prior distribution will be omitted from the model description. The model description can be manually edited.


Refer to SimPol in the Literature for real data examples of how model comparison has been used in inference.

Rendering elements onto the screen (for example polymerase motion, plots, simulated sequences) can be computationally costly and can slow down the whole program. If you are running a complex simulation on a long sequence over a large number of trials, you may notice the simulation begin to slow down eventually. Maximising the performance of this program requires you as the user to regulate when rendering does or does not occur. If the program starts to slow down, you should hide components of the interface with its respective minimise − button. If a component is minimised then rendering will not occur. The polymerase animation can be minimised by setting the speed to hidden. This will happen automatically when simulating large sequences.

The simulation component of SimPol is written in C++ and compiled with WebAssembly. You will be notified if your web browser fails to instantiate the WebAssembly module. In this case, the simulation will instead be run using the JavaScript module which will have a significant impact on performance.

In January 2018 we tested SimPol on different operating systems using recent versions of different web browsers, using both JavaScript (JS) and C++. Each setting was tested on 1000 simulations of the 4029nt E.coli ropB gene (accession number KP670789). Compiling the C++ code with WebAssembly increases speed performance approximately 10x over the same algorithm written in JavaScript for the fastest JavaScript environments. For maximum performance we recommend using SimPol on Firefox, Chrome, Safari, or through the command line. Benchmark session: Open | View XML

OS	Device description	Program (SimPol language)	Time (min)
Ubuntu 7	Intel i5 3.30GHz	Command prompt (C++)
		Node.js 8.6 (JS)	3:02
		Chrome 61 (C++)	0:15
		Chrome 61 (JS)	3:22
		Firefox 58 (C++)	0:14
		Firefox 58 (JS)	5:14
Windows 10	AMD Ryzen 5 1600 3.20 GHz	Command prompt (C++)
		Cygwin (C++)	0:37
		Node.js 8.7 (JS)	2:49
		Chrome 63 (C++)	0:15
		Chrome 63 (JS)	2:37
		Firefox 58 (C++)	0:16
		Firefox 58 (JS)	4:37
		Microsoft Edge 40 (C++)
		Microsoft Edge 40 (JS)	15:33
OS X 10.13	Intel i5 3.0GHz	Command prompt (C++)	0:14
		Safari 10.3 (C++)	0:13
		Safari 10.3 (JS)	2:38
iOS 10.0.3	iPhone 7	Chrome 61 (JS)
		Safari 10.3 (JS)
Android 7.0	Samsung Galaxy S6 Edge+	Chrome 64 (C++)	0:52
		Chrome 64 (JS)	8:03
		Firefox 58 (C++)	0:46
		Firefox 58 (JS)	11:23

The following instructions will work for Linux/Mac.

In this project (publication under review), we inferred kinetic parameters for three RNA polymerases: E. coli RNAP, S. cerevisiae pol II, and Bacteriophage T7 pol. The posterior distributions for these three experiments are available below. Posterior distributions may be viewed directly in the web browser or downloaded. Downloaded .log files can also be viewed with Tracer.

Enzyme	View .xml	Model	Open	Download .log
E. coli RNAP	Ecoli_RNAP_session.xml	11	Open	Download
		12	Open	Download
S. cerevisiae pol II	Yeast_polII_session.xml	11	Open	Download
		12	Open	Download
Bacteriophage T7 pol	T7_pol_session.xml	5	Open	Download