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).

Variables

Variable	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_{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	\(0.5\)
\(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	\(10^{-5}\)	[8][12]
\(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

We have used a small observation box to represent the whole cell. 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} \]

The diffusion coefficient for a single monomer inside the cell is \(D_{cell}=4.6( \mu m)^{2}/ms \) ^[14] and for the outside of the cell, \(D_{outside}=5( \mu m)^{2}/ms \), with a continuous intermediate for the membrane of the cell.

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_+~10^5 mol^{-1} s^{-1} \) and the fragmentation rate to be of order \(k_-~10^{-8} s^{-1} \), the rate of elongation is much higher than that of fragmentation as \( \frac{k_+}{k_-}-=10^{13} mol^{-1} \).

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 angle \(\alpha\) is the angle between the tangent of the first monomer of the chain and the tangent of the monomer at distance \(L \) away on the chain. The persistence length of the chains determines the angles at which the free monomers can bind. The persistence length 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]

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:

\[ <\cos (\alpha)> = \exp \big(- \frac{L}{\lambda } \big) = \exp \big(- \frac{jl}{\lambda} \big) \] \[ 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.

Chains detaching from the membrane

Each chain 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.

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.