Team:Kent/Modeling

Modeling

Modeling, i.e. the mathematical description of a physical system, is an important aspect of the scientific method as it allows us to quantitatively understand the system and to make predictions. Specifically, a model helps understanding how changes in any parameter affect the system; this is especially important when the relevance of some parameters are unknown. We chose to develop a simulation model of the system studied in our project, which consists of a self-assembling biological structure that we thought would be exciting and informative to quantify and visualize in an interactive way. We made the computer code of the model publicly available, so that other teams can play around and build upon it.

We used the Monte Carlo method to simulate the diffusion ^[3] of monomers inside an E.coli cell, their transit through the cell membrane and the production of amyloid nanowires outside the cell. A typical E.coli cell has a length, \(l=2 \mu m \) and a diameter, \(d=1 \mu m \). We take a small observation cube, which is a portion of both the inside of the cell and the bulk outside the cell; we can extrapolate this to describe the whole system. Simulating just a portion has the advantage of requiring less computational power.

The monomers initially start at the bottom of the observation volume and are allowed to stochastically diffuse inside the cell and go through the cell membrane. The binding site seeds are located on top of the cells membrane and when a monomer gets close enough to the binding site it may bind and form a link in the chain. Over time the chains can grow to lengths in the range of \(60nm \) to \(100 \mu m \)^[1][5][10] and have diameter \(20nm\) ^[13][14], these chains are not necessarily straight and persistence length parameter dictates the angles at which particles can bind. The monomer can only interact with the monomer attached to the chain, i.e. there is no branching^[13].

Periodic and reflective boundary conditions are applied in the model. When a particle leaves through the side of the observation volume, we can assume that another particle enters through the other side. When a particle reaches the bottom of the observation volume we can assume that the particle is reflected. This allows us to reproduce behaviour similar to the bulk.

The number of monomers in the system is not necessarily constant; monomers can be created and degraded. However, Hall ^[8] (2003) argued with the two extreme cases that either; the monomer are generated and degraded on a time scale much slower than amyloid growth so the number of particles are constant; or that the amyloid growth occurs on a time scale much slower than monomer degradation so we can refer to the free concentration of the monomer as constant.

The Monte Carlo simulation was implemented in Matlab and visualized using Visual Molecular Dynamics (VMD). Our software, including our Matlab code can be found here: Team Kent Software

Parameters

Parameters	Description	Value	Reference
\(L_x\)	Length of the observation box volume	\(0.4 \mu m\)
\(L_y\)	Depth of the observation box volume	\(0.4 \mu m\)
\(L_z\)	Height of the cell in observation volume	\(0.5 \mu m\)
\(N_b\)	Number of binding sites	\(5\)
\(N(0)\)	Initial number of particles	2000
\(\Delta t\)	Length of each timestep	\(\frac{1}{2000} s\)
\(\lambda \)	Persistence Length	\(23 \mu m\)	[7][11]
\(r\)	Monomer radius	\( 10nm\)	[13] [14]
\(r_b \)	Binding radius	\(20nm \)	[13] [14]
\(D_{cell} \)	Diffusion coefficient of the monomer inside the cell	\(4.6 \mu m^2 /s \)	[14]
\(D_{outside} \)	Diffusion coefficient of the monomer inside the cell	\(30 \mu m^2 /s \)	[14]
\(D_{chain} \)	Diffusion coefficient of the chain after the chain has detached	\(0.5 \mu m^2 /s \)
\(K_{+} \)	Reaction rate of elongation	\(10^{5} mol^{-1}s^{-1}\)	[8][12]
\(K_{-} \)	Reaction rate of fragmentation	\(10^{-8} s^{-1}\)	[8][12]
\(P(bind) \)	Probability that a monomer, within range, will attach to the end of a chain	Variable
\(P(enter) \)	Probability that a monomer will enter the cell from the outside	\(0.5\)
\(P(leave) \)	Probability that a monomer will leave the cell from the inside	\(0.5\)
\(P(spawn) \)	Probability that a monomer will be created at the end of each timestep	\(0\)	[8]
\(P(decay) \)	Probability that any given monomer will decay at each timestep	\(0\)	[8]
\(P(es) \)	Probability that the end monomer on a chain will detach	Variable
\(P(anchor detach) \)	Probability that the chain will detach from the membrane	\(10^{-6}\)
\(P(is) \)	Probability that a chain of \(j\) monomers will fragment internally	\( (j-1)*P(es) \)	[8][12]

Observation box volume and area

Figure 1: A simple 2D schematic of the cell, bulk and observation box

We have used a small observation box to represent the whole cell, (see figure 1). It is important to know what proportion of the cell is represented so that the results can be extrapolated properly. Only particles represented in the observation box can interact with the system.

The typical length of an E.coli cell is \(L \approx 2 \mu m\) and the typical diameter, \(d \approx 1 \mu m\). If we consider the cell to be composed of a cylinder of length \(l=L-d\) and two hemispheres at each end then the volume of the sphere is: \[V_{cell} = V_{cylinder} + V_{sphere}\] \[V_{cell} = \pi \big(\frac{d}{2}\big)^2 (L-d)+ \frac{4}{3} \pi \big(\frac{d}{2}\big)^3 = 1.309 (\mu m)^3\] The volume of the cell inside the observation box is \[V_{ob} = L_x L_y L_z = 0.08(\mu m)^3\] From this we can work out the proportion of the volume of cell inside the observation box \[\frac{V_{ob}}{V_{cell} } = 6.112 \% \] Likewise for the surface area of the cell \[A_{cell} = A_{cylinder} + A_{sphere}\] \[A_{cell} = \pi d(L-d) + 4 \pi \big(\frac{d}{2}\big)^2 = 6.28 (\mu m)^2 \] \[A_{ob} = L_x L_y = 0.16 (\mu m)^2 \] The proportion of the cell's surface area inside the observation box is:

\[\frac{A_{ob}}{A_{cell}} = 2.54\% \]

Stochastic Brownian motion

The displacement each time step in the \(x\),\(y\) and \(z\) directions is given by a simple form of the Langevin equation. This model assumes that there are no deterministic forces acting on the monomers, i.e. pure diffusion ^[3][4].

\[ \Delta x = \xi \sqrt{2D \Delta t} \]

Kawai-Noma ^[14] found the diffusion coefficient for a single Sup35 monomer inside the cell is \(D_{cell}=4.6( \mu m)^{2}/s \). We rounded up the theoretical diffusion coefficient that Kawai-Noma ^[14] calculated from the molecular mass to get the diffusion coefficient for the outside of the cell, \(D_{outside}=30( \mu m)^{2}/s \).

Boundary conditions

We make the assumption that if a monomer leaves through the side of the observation volume that another particle will enter through the other side. We also assume that if a monomer hits the bottom of the observation volume that it will be reflected. In the following equations, \(x' \), \(y' \) and \(z' \) denote the respective \(x\),\(y\) and \(z\) coordinates of the particles after the application of boundary conditions:

\[ x' = \begin{cases} x+L_x & x < 0 \\ x & 0 < x < L_x \\ x - L_x & L_x < x \end{cases} \] \[ y' = \begin{cases} y+L_y & y < 0 \\ y & 0 < y < L_y \\ y - L_y & L_y < y \end{cases} \] \[z' = -z \]

If a particle on the inside moves to a position that has a height greater than the height of the cell then it has a probability \(P(leave) \) of successfully leaving the cell. Likewise if a particle on the outside of the cell moves to a position that has a height less than the height of the cell then it has a probability \(P(enter) \)of successfully entering the cell. In both of these scenarios, if the particle fails to pass through the membrane then they are relocated in the \(z \) position by:

\[z' = 2 L_z - z \]

Length of the chain

There are three processes which may affect the length of a chain, nucleation, elongation; and fragmentation. For each chemical equation \(A \) represents a monomer, \(B \) represents a binding site; \(AB \) represents a monomer-binding site complex and \(A_j B \) denotes a chain of \(j \) monomers attached to a binding site. For the reaction rates, a least squares fit from both Hall^[8](2004) and Xue^[12](2013), revealed the elongation rate to be of order \(k_+ \approx 10^5 mol^{-1} s^{-1} \) and the fragmentation rate to be of order \(k_- \approx 10^{-8} s^{-1} \).

Each chain also has a probability of detaching from the membrane, once this event occurs the whole chain, including the original binding site, is free to diffuse in the bulk. This is similar to the diffusion of a fragment, wherein the central monomer of the chain will diffuse and the other particles of the chain will remain relative to it.

Nucleation

\[ A + B \overset{k}{\rightarrow} AB\]

For simplicity we have made no distinction between nucleation and elongation.

Elongation

\[ A_{j}B \overset{k_{+}}{\rightarrow} A_{j + 1}B \]

If a free monomer is close enough, a distance less than \(r_b\), to a monomer-binding site, it may bind. However, if the monomer is too close to the monomer-binding site, within the radius of the monomer \(r\), the monomer-binding site will repel the free monomer.

The persistence length of the chains determines the angles at which the free monomers can bind. The persistence length, \(\lambda\), of long amyloid fibrils with near 100% \(\beta \) -sheet content is about \(23 \mu m\), roughly 40 times higher than the persistence length of short worm-like chains. ^[7][11] The angle \(\alpha\) is the angle between the tangent of the first monomer of the chain (at position \(0\)) and the tangent of the \(j\)th monomer at distance \(L \) away from the first monomer. The expectation value of the cosine of the angle falls off exponentially as the distance from the first polymer increases. \[ <\cos (\alpha)> = \exp \big(- \frac{L}{\lambda } \big) = \exp \big(- \frac{jl}{\lambda} \big) \] Where \(l\) is the length of a single monomer.

Fragmentation

\[ A_{j}B \overset{k_{-}}{\rightarrow} A_{i}B + A_{j-i} \]

In end scission the last monomer on a chain breaks away from the rest of the chain. In internal scission a group of \(i \) particles break away from a chain of j polymers ( \( i < j \) ). The probability of an internal scission event occurring increases linearly with the length, \(P(is)=(j-1)*P (es) \) of the chain due to a chain of more polymers has more links that could spontaneously break. ^[8]

When a fragment breaks off of a chain, the fragment can diffuse freely in the bulk. The position of the middle particle in the fragment is determined using the Langevin equation and the rest of the particles in the fragment are positioned relative to the middle particle using spherical coordinates.

We calculate the relation of all the particles in the fragment in spherical coordinates compared to the center monomer of the fragment, where:

\[ d_x = r \sin(\theta) \cos(\phi) \] \[ d_y = r \sin(\theta) \sin (\phi) \] \[ d_z = r \cos (\phi) \]

We then determine the position of the central particle and then find the position of the other particles in the chain relative to it, using

\[ x = x_c + r \sin (\theta + \Delta \theta ) \cos (\phi + \Delta \phi \Delta t) \] \[ y = y_c + r \sin (\theta + \Delta \theta) \sin (\phi + \Delta \phi \Delta t) \] \[ z = z_c + r \cos (\theta + \Delta \theta ) \]

Where \(x_c \), \(y_c \) and \(z_c \) are the coordinates of the central particle, \(x\), \(y\) and \(z\) are coordinates of the other particles in the fragment. \( \Delta \theta \) and \( \Delta \phi \) are constants and \(n \) is the number of timesteps since the chain detached from the membrane.

Results

After creating the simulation, we changed some of the variables and recorded how this affected the system. Graphs 1 and 2 show the effect of changing the probability that a monomer will bind when within range of a binding site. In graph 1 we see that the binded particle to total particle ratio increases over time, the gradient at the beginning is steep but begins to lessen over time. The simulations with a higher \(P(bind) \) have a much steeper gradient than those with a lesser value of \(P(bind) \).

In graph 2 we see that for each simulation, the mean length of the fibrils increases over time with a steep gradient that lessens over time. At \(1500 ms\) the mean length of the chains for the \(P(bind) = 1\) simulation is over twice that of the mean length of the chains corresponding to the \(P(bind) = 0.25 \) simulation.

In graph 3 all the simulations show an initial rapid increase of mean chain length. The \(P(es) = 10^{-6} \) and \(P(es) = 10^{-5} \) continue this trend with a gradual lessening of the gradient for the duration of the simulation, however the mean length of the chain in the \(P(es) = 10^{-4} \) and \(P(es) = 10^{-3} \) simulations start to fluctuate. The fluctuation range for \(P(es) = 10^{-4} \) is much greater than that for \(P(es) = 10^{-3} \). The \(P(es) = 10^{-6} \) and \(P(es) = 10^{-5} \) mean chain length values are close to each other until at around \(1250ms \) the mean chain length for the \(P(es) = 10^{-5} \) simulation is cut by about 20 polymers. The mean length of the chain in the \(P(es) = 10^{-4} \) simulation is consistently higher than the mean length of the chain in the \(P(es) = 10^{-3} \) simulation.

Graph 1: The ratio of binded monomers to free monomers over time for different values of \(P(bind) \). \(P(es) \) is set to \(10^{-5} \)

Graph 2: The mean length of the chains over time for different values of \(P(bind) \). \(P(es)\) is set to \(0\).

Graph 3: The mean length of the chains over time for different values of \(P(es)\). \(P(bind) \) is set to \(0.5 \).

Conclusions

In conclusion, we have made a simulation of the self-assembling structure studied in our project using Monte Carlo Brownian dynamics. To ensure maximum accuracy of the model, real world values were added to the system. We have made the computer code for our simulation publicly available and have included a graphical user interface to encourage other teams to play around and build upon it. We then explained some of the key features of our model such as elongation and fragmentation. Finally, we changed some parameters of our simulation and recorded how these changes affected the system.

References

[1] Sivanathan, V., & Hochschild, A. (2012). Generating extracellular amyloid aggregates using E. coli cells. Genes & development, 26(23), 2659-2667.

[2] Elowitz, M. B., Surette, M. G., Wolf, P. E., Stock, J. B., & Leibler, S. (1999). Protein Mobility in the Cytoplasm of Escherichia coli. Journal of bacteriology,181(1), 197-203.

[3] Philipse, A. P. (2011). Notes on Brownian Motion. Utrecht University, Debye Institute, Van’t Hoff Laboratory.

[4] Barlett, V. R., Hoyuelos, M., & Mártin, H. O. (2013). Monte Carlo simulation with fixed steplength for diffusion processes in nonhomogeneous media. Journal of Computational Physics, 239, 51-56.

[5] Xue, W. F., Homans, S. W., & Radford, S. E. (2008). Systematic analysis of nucleation-dependent polymerization reveals new insights into the mechanism of amyloid self-assembly. Proceedings of the National Academy of Sciences,105(26), 8926-8931.

[6] Knowles, T. P., Waudby, C. A., Devlin, G. L., Cohen, S. I., Aguzzi, A., Vendruscolo, M., ... & Dobson, C. M. (2009). An analytical solution to the kinetics of breakable filament assembly. Science, 326(5959), 1533-1537.

[7] Smith, J. F., Knowles, T. P., Dobson, C. M., MacPhee, C. E., & Welland, M. E. (2006). Characterization of the nanoscale properties of individual amyloid fibrils.Proceedings of the National Academy of Sciences, 103(43), 15806-15811.

[8] Hall, D., & Edskes, H. (2004). Silent prions lying in wait: a two-hit model of prion/amyloid formation and infection. Journal of molecular biology, 336(3), 775-786.

[9] Lomakin, A., Chung, D. S., Benedek, G. B., Kirschner, D. A., & Teplow, D. B. (1996). On the nucleation and growth of amyloid beta-protein fibrils: detection of nuclei and quantitation of rate constants. Proceedings of the National Academy of Sciences, 93(3), 1125-1129.

[10] Scheibel, T., Parthasarathy, R., Sawicki, G., Lin, X. M., Jaeger, H., & Lindquist, S. L. (2003). Conducting nanowires built by controlled self-assembly of amyloid fibers and selective metal deposition. Proceedings of the National Academy of Sciences, 100(8), 4527-4532.

[11] vandenAkker, C. C., Engel, M. F., Velikov, K. P., Bonn, M., & Koenderink, G. H. (2011). Morphology and persistence length of amyloid fibrils are correlated to peptide molecular structure. Journal of the American Chemical Society,133(45), 18030-18033.

[12] Xue, W. F., & Radford, S. E. (2013). An imaging and systems modeling approach to fibril breakage enables prediction of amyloid behavior. Biophysical journal, 105(12), 2811-2819.

[13] Glover, J. R., Kowal, A. S., Schirmer, E. C., Patino, M. M., Liu, J. J., & Lindquist, S. (1997). Self-seeded fibers formed by Sup35, the protein determinant of [PSI+], a heritable prion-like factor of S. cerevisiae. Cell, 89(5), 811-819.

[14] Kawai-Noma, S., Pack, C. G., Kojidani, T., Asakawa, H., Hiraoka, Y., Kinjo, M., ... & Hirata, A. (2010). In vivo evidence for the fibrillar structures of Sup35 prions in yeast cells. The Journal of cell biology, 190(2), 223-231.