ReaDDy - A particle-based
reaction-diffusion simulator

Practical sessions

Here we’ll gather the material and tasks for the daily sessions.

Solutions to tasks we have already done, will be uploaded here

Monday - ReaDDy intro

Task 0) installation

Get miniconda and install it

wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh

Check your quota -s. If your home directory is limited in space, you might want to install it under storage, i.e. when prompted for install location, choose either /home/mi/user/miniconda3 or /storage/mi/user/miniconda3 (replace user with your username). Your ~/.bashrc should contain the line

. /storage/mi/user/miniconda3/etc/profile.d/conda.sh

to enable the conda command. If you try out which conda, you should see its function definition, and you’re good to go.

Create a conda environment workshop

conda create -n workshop

and activate it.

conda activate workshop

Add conda-forge as default channel

conda config --add channels conda-forge --env

and install the latest readdy in it

(workshop) $ conda install -c readdy/label/dev readdy

Check if it worked, start a python interpreter and do

import readdy
print(readdy.__version__)

If this does not return an error, you are ready to readdy (almost). Make sure you also have vmd installed. It should be installed on the universities machines.

Task 1) Python, ipython notebook, numpy, matplotlib crash course

Follow along the ipython-intro. You should be able to reproduce the presented usage in your own ipython notebook. The notebook covers the following:

  • usage of an ipython notebook
  • filesystem operations with os
  • save and load data objects with pickle
  • numerical operations with numpy
  • plotting with matplotlib.pyplot

Task 2) ReaDDy intro I: Particles, diffusion and potentials

Follow along the readdy-intro. You should be able to reproduce the presented usage in your own ipython notebook. The notebook covers the following:

  • principle workflow of readdy: system and simulation
  • adding particle species to system
  • adding potentials to system
  • spatial layout of simulation box and box potentials
  • adding particle instances at given positions to simulation
  • convert trajectory output to VMD viewable format

Task 3) ReaDDy intro II: Reactions, observables and checkpoints

Follow along the readdy-observables The notebook covers the following:

  • adding reaction to system
  • adding observable to simulation
  • reading back the observable from trajectory file
  • using checkpoints to continue a simulation

Task 4) Crowded mixture, MSD

The time dependent mean squared displacement (MSD) is defined as

$$ \langle(\mathbf{x}_t-\mathbf{x}_0)^2\rangle_N $$

where $\mathbf{x}_t$ is the position of a particle at time $t$, and $\mathbf{x}_0$ is its initial position. The difference of these is then squared, yielding a squared travelled distance. This quantity is then averaged over all particles.

The task is to set up a crowded mixture of particles and measure the MSD. There is one species A with diffusion coefficient 0.1. The simulation box shall be $20\times20\times20$ and non-periodic, i.e. define a box potential that keeps particles inside, e.g. with an extent $19\times19\times19$ and appropriate origin.

Define a harmonic repulsion between A particles with force constant 100 and interaction distance 2.

Add 1000 A particles to the simulation, uniformly distributed within the box potential.

Observe the positions of particles with a stride of 1.

Run the simulation with a time step size of 0.01 for 20000 steps.

In the post-processing you have to calculate the MSD from the observed positions. Implement the following steps:

  1. Convert positions list into numpy array $T\times N\times3$
  2. From every particles position subtract the initial position
  3. Square this
  4. Average over particles np.mean(...)

Since the positions observable returns a list of list, it might come in handy to convert this to a numpy array of shape $T\times N\times3$ (step 1). You may use the following brute force method to do this

T = len(positions)
N = len(positions[0])
pos = np.zeros(shape=(T, N, 3))
for t in range(T):
    for n in range(N):
        pos[t, n, 0] = positions[t][n][0]
        pos[t, n, 1] = positions[t][n][1]
        pos[t, n, 2] = positions[t][n][2]

You shall implement steps 2.-4. by yourself.

4a) Finally plot the MSD as a function of time in a log-log plot, let’s call this the crowded result. (Hint: You may find that there is an initial jump in the squared displacements. Equilibrating the particle positions before starting the measurement may help you there. Make use of checkpointing to use an already equilibrated state)

4b) Repeat the whole procedure, but do not register the harmonic repulsion between particles, this shall be the free result. Compare the MSD to the previous result, possibly in the same plot.

4c) Additionally plot the analytical solution for freely diffusing particles

$$ \langle(\mathbf{x}_t-\mathbf{x}_0)^2\rangle_N = 6 D t $$

From your resulting plot identify “finite size saturation”, “subdiffusion”, “reduced normal diffusion”, “ballistic/normal diffusion”

Task 5) Crowded mixture, RDF

The radial distribution function (RDF) $g(r)$ describes how likely it is to find two particles at distance $r$. Compare the RDF of harmonically repelling particles and the RDF of non-repelling particles.

Therefore set up the same system as in Task 4) but this time the system shall be periodic and there is no box potential.

You may want to equilibrate the initial particle positions, use checkpointing.

Instead of observing the particle positions, observe the radial distribution function

simulation.observe.rdf(stride=1000, bin_borders=np.arange(0.,10.,0.2),
                       types_count_from="A", types_count_to="A",
                       particle_to_density=n_particles/system.box_volume)

In the post-processing, you shall use

_, bin_centers, distribution = traj.read_observable_rdf()

to obtain the observable. The distribution contains multiple $g(r)$, one for each time it was recorded. Average them over time.

5a) Plot the RDF as a function of the distance (i.e. meandistribution as function of bin_centers)

5b) Perform the same procedure but for non-interacting particles, compare with the previous result, preferably in the same plot. What are your observations?

Tuesday

This session will cover a well studied reaction diffusion system, the Lotka Volterra system, that describes a predator prey dynamics. We will investigate how parameters in these system affect spatiotemporal correlations.

Task 1) Diffusion influenced reaction

Consider the reaction

$$ A + A\to B\quad\text{with rate: }\lambda=0.1 $$

The reaction distance should be $R=1$.

Set up a periodic box with dimensions $20\times20\times20$

Add 1000 A particles with $D=1$. The diffusion constant of B should be $D=1$. They should be uniformly distributed in the box, use

n_particles = 1000
init_pos = np.random.uniform(size=(n_particles,3)) * extent + origin
simulation.add_particles("A", init_pos)

Observe the number of particles with a stride of 1. Additionally you can print the current number of particles using the callback functionality

simulation.observe.number_of_particles(
  1, ["A","B"],
  callback=lambda x: print("A {}, B {}".format(x[0],x[1]))
)

Simulate the system for 10000 steps with time step size 0.01.

1a) Obtain the number of particles and plot them as a function of time.

The analytic solution of the concentration of particles subject to the reaction above is obtained by solving the ODE

$$ \frac{\mathrm{d}a}{\mathrm{d}t} = -k a^2\quad a(0)=a_0 $$

which yields

$$ a(t) = \frac{1}{a_0^{-1}+kt} $$

1b) Fit the function $a(t)$ to your counts data, to obtain the parameter $k$. For $a_0$ use n_particles. Make use of scipy.optimize.curve_fit. Can the function $a(t)$ describe the data well?

1c) Repeat the procedures a) and b) but change the diffusion coefficient to $D=0.01$ and change the microscopic reaction rate to $\lambda=1$. What happened to your fit?

Task 2) well mixed predatory-prey

Simulate a Lotka-Volterra/predator-prey system for the given parameters and determine the period of oscillation.

There are two particle species A (prey) and B (predator) with the same diffusion coefficient $D=0.7$. Both particle species shall be confined to a quadratic 2D plane with an edge length of roughly 50, and non periodic boundaries, i.e. the 2D plane must be fully contained in the simulation box. Choose a very small thickness for the plane, e.g. 0.02, and a force constant of 50.

Particles of the same species interact via harmonic repulsion, with a force constant of 50, and an interaction distance of 2.

There are three reactions for the Lotka Volterra system: birth, eat, and decay.

$$ \mathrm{birth}: \mathrm{A}\to \mathrm{A} +\mathrm{A}\quad\mathrm{with }~ \lambda=4\times 10^{-2}, R=2 $$
$$ \mathrm{eat}: \mathrm{A}+\mathrm{B}\to \mathrm{B} +\mathrm{B}\quad\mathrm{with }~ \lambda=10^{-2}, R=2 $$
$$ \mathrm{decay}: \mathrm{B}\to \emptyset\quad\mathrm{with }~ \lambda=3\times 10^{-2} $$

Initialize a system by randomly positioning 250 particles of each species on the 2D plane. Run the simulation for 100,000 integration steps with a step size of 0.01. Observe the number of particles and the trajectory, with a stride of 100.

From the observable, plot the number of particles for the two species as a function of time in the same plot. Make sure to label the lines accordingly. Additionally plot the phase space trajectory, i.e. plot the number of predator particles against the number of prey particles.

Question: What is the period of oscillation?

Task 3) diffusion influenced predator-prey

We simulate the same system as in Task 2) but now introduce a scaling parameter $\ell$. This scaling parameter is used to control the ratio of reaction rates $\lambda$ and the rate of diffusion $D$.

For a given parameter $\ell$

  • change all reaction rate constants $\tilde\lambda=\lambda\times\sqrt{\ell}$.
  • change all diffusion coefficients $\tilde D= D/\sqrt{\ell}$.

where the tilde ~ denotes the actually used value in the simulation and the value without tilde is the one from Task 1).

This means that the ratio

$$ \frac{\tilde\lambda}{\tilde D} = \ell\,\frac{\lambda}{D}=\ell\times\mathrm{const} $$

is controlled by $\ell$.

Perform the same simulation as in Task 1) but for different scaling parameters $\ell$, each time saving the resulting plots (time series and phase plot) to a file, use plt.savefig(...). Run each simulation for $n$ number of integration steps, where

$$ n(\ell) = 10^4 + \frac{10^5}{\sqrt{\ell}} $$

Vary the scaling parameter from 1 to 400.

Question: For which $\ell$ do you see a qualitative change of behavior in the system, both from looking at the plots and the VMD visualization? Which cases are reaction-limited and which cases are diffusion-limited?

Task 4)

Starting from the parameters of a diffusion influenced system from Task 3), set up a simulation, where prey particles form a traveling wavefront, closely followed by predators. Therefore you might want to change the simulation box and box potential parameters to get a 2D rectangular plane.

Feel free to adjust all parameters. You can experiment with other spatial patterns as well, have a look at this paper.

Wednesday

This session will cover macromolecules and their dynamics. In particular we want to model short RNA chains, represented by linear chains of beads.

Assume $N$ beads of particles at positions $\mathbf{x}_i$, the $i$th bead is connected with the $i+1$th bead by a spring which has a fixed length $l$. Thus the whole chain of particles is linearly connected. The vector $\mathbf{r}_i=\mathbf{x}_{i+1}-\mathbf{x}_i$ is the segment that connnects adjacent beads.

One can measure how strongly the $i$th and the $j$th segment correlate by considering the scalar product $\mathbf{r}_i\cdot\mathbf{r}_j$.

Since this value alone is quite meaningless, one can consider its average over a whole ensemble of segments. This average can be taken over all segments, i.e. $\forall i,j$ in the linear chain, and also over many times if the linear chain evolves over time. For realistic polymers one typically observes that the $i$th segment strongly correlates with the $j=i+1$th segment. However for $j\gg i$ the correlation vanishes. Phenomenologically this is understood by an exponential decay of the correlation

$$ \langle \mathbf{r}_i\cdot\mathbf{r}_j \rangle=l^2\exp\left(-|j-i|\,l/l_p\right) $$

where we have defined the persistence length $l_p$. The value of $l_p$ is determined by the interaction of the beads. For example structured RNA molecules typically show a persistence length of 72nm. Today we will instead model unstructed RNA molecules, which show a persistence length of roughly 2nm.

We will consider two models for the polymer, namely

  • the freely jointed chain (FJC), and
  • the freely rotating chain (FRC).

In the FJC, the beads are connected by segments of fixed length $l=0.48$. Other than that there is no interaction.

In the FRC, the beads are also connected by segments of fixed length $l=0.48$. Additionally the angle between neighbouring segments is fixed to $\theta=35^\circ$.

Task 1) Equilibrate polymers

In the first task you will

  • a) equilibrate a freely jointed chain (FJC) of $N=50$ beads. Equilibration is ensured by measuring the radius of gyration of a polymer over time.
  • b) Once the polymer is equilibrated, you will measure the persistence length $l_p$.
  • c) Repeat a) and b) for the freely rotating chain (FRC)

We will only need one particle species monomer with diffusion constant $D=0.1$. Note that this will not be a normal species but will be a topology species:

system.add_topology_species("monomer", 0.1)

The simulation box shall have box_size= [102.,102.,102.], and be non-periodic. This means that there has to be a box potential registered for the type monomer with force constant 50, origin=np.array([-50,-50,-50]) and extent=np.array([100,100,100]).

In order to build a polymer we use topologies. At first we need to register a type of toplogy

system.topologies.add_type("polymer")

The monomers in a polymer must be held together by harmonic bonds, defined by pairs of particle types

system.topologies.configure_harmonic_bond(
    "monomer", "monomer", force_constant=50., length=0.48)

Only in the case of a FRC we also specify an angular potential, that is defined for a triplet of particle types

system.topologies.configure_harmonic_angle(
    "monomer", "monomer", "monomer", force_constant=20., 
    equilibrium_angle=(180.-35.)/180.*np.pi)

where an equilibrium_angle is given in radians (Note that the equilibrium angle here is not the same as the $\theta$ as defined above, thus the conversion by 180 degrees).

The next step is creating the simulation. Then we specify the observables, the trajectory and the particle positions

simulation.record_trajectory(stride=10000)
simulation.observe.particle_positions(stride=1000)

We also want to make use of checkpointing to continue simulation for an already equilibrated polymer. If there are no checkpoints, we want to create new positions for the polymer. The new positions represent a random walk in three dimensions with fixed step length bond_length.

if os.path.exists(checkpoint_dir):
    # load checkpoint
    simulation.load_particles_from_latest_checkpoint(checkpoint_dir)
else:
    # new positions
    init_pos = [np.zeros(3)]
    for i in range(1, chain_length):
        displacement = np.random.normal(size=(3))
        displacement /= np.sqrt(np.sum(displacement * displacement))
        displacement *= bond_length
        init_pos.append(init_pos[i - 1] + displacement)
    init_pos = np.array(init_pos)

    # subtract center of mass
    init_pos -= np.mean(init_pos, axis=0)

    # add all particles for the topology at once
    top = simulation.add_topology("polymer", len(init_pos) * ["monomer"], init_pos)

    # set up the linear connectivity
    for i in range(1, chain_length):
        top.get_graph().add_edge(i - 1, i)

# this also creates the directory
simulation.make_checkpoints(n_steps // 100, 
  output_directory=checkpoint_dir, max_n_saves=10)

Tip: Keep two separate checkpoint directories output files and for the FJC and the FRC model. This means you may want to have the following defined in the beginning of your notebook

chain_type = "fjc" # fjc or frc
out_dir = "/some/place/on/your/drive"
out_file = os.path.join(out_dir, f"polymer_{chain_type}.h5")
checkpoint_dir = os.path.join(out_dir, f"ckpts_{chain_type}")

Now that we have defined the simulation object we can run the simulation

if os.path.exists(simulation.output_file):
    os.remove(simulation.output_file)
simulation.run(n_steps, dt)

and also observe the output

traj = readdy.Trajectory(out_file)
traj.convert_to_xyz(
  particle_radii={"monomer": bond_length / 2.},
  draw_box=True)

1a) The radius of gyration is a measure of how ‘extended’ in space a polymer is. To calculate it, we must have observed the particle positions. As a first step we convert the readdy output to a numpy array

times, positions = traj.read_observable_particle_positions()

# convert to numpy array
T = len(positions)
N = len(positions[0])
pos = np.zeros(shape=(T, N, 3))
for t in range(T):
    for n in range(N):
        pos[t, n, 0] = positions[t][n][0]
        pos[t, n, 1] = positions[t][n][1]
        pos[t, n, 2] = positions[t][n][2]

Then from the pos array you may use the following to calculate the radius of gyration (the assertion statements may help you understand how the arrays are shaped)

# calculate radius of gyration
mean_pos = np.mean(pos, axis=1)

difference = np.zeros_like(pos)
for i in range(n_particles):
    difference[:,i] = pos[:,i] - mean_pos

assert difference.shape == (T,N,3)

# square and sum over coordinates (axis=2)
squared_radius_g = np.sum(difference * difference, axis=2)

assert squared_radius_g.shape == (T,N)

# average over particles (axis=1)
squared_radius_g = np.mean(squared_radius_g, axis=1)

radius_g = np.sqrt(squared_radius_g)

assert radius_g.shape == times.shape == (T,)

Plot the radius of gyration as a function of time, is it equilibrated? If not, simulate for a longer time.

1b) The mean correlation of segments $\langle \mathbf{r}_i\cdot\mathbf{r}_j \rangle$ shall be calculated from the pos array. You will average over all pairs of the linear chain and also over all times.

Use the following snippet to calculate the segments vector

assert pos.shape == (T, N, 3)
# calculate segments
segments = pos[:, 1:, :] - pos[:, :-1, :]

The correlation between $i$ and $j$ shall be measured as a function of their separation $s=\mid j-i\mid$, which is a value between 0 and $N-1$, e.g.

n_beads = pos.shape[1]
separations = np.arange(0, n_beads - 1, 1)
corrs = None # your task

Now for every separation $s$, calculate the average correlation, averaged over all pairs $(i,j)$ that lead to $s=\mid j-i\mid$.

Hints:

  • The calculation may involve a double loop over all segments

    for i in range(n_beads-1):
        for j in range(i, n_beads-1):
            # something
    
  • The scalar product of the $i$th and the $j$th bead for all times is

    np.sum(segments[:, i, :] * segments[:, j, :], axis=1)
    

    where the summation over axis=1 is over the x,y,z coordinates

Then plot the mean correlation as a function of the separation. What do you observe?

To determine the persistence length $l_p$, fit a function of the form

$$ f(s)=c_1e^{-sc_2} $$

to the data using scipy.optimize.curve_fit. From the constant $c_2$ you should be able to obtain the persistence length. What is its value?

1c) Repeat a) and b) for the FRC. Note the value for $l_p$ for both models.

Given the values in the introduction text, which of the models (FJC, FRC) is better suited to model unstructured RNA?

Task 2) Identify given structures

Obtain the two data files a.npy and b.npy and save them to a directory of your liking. Each of them contains 100 positions of beads, i.e. a and b are two polymer configurations. You can load them as follows

a = np.load("a.npy")
b = np.load("b.npy")
assert a.shape == (100, 3)
assert b.shape == (100, 3)

Your task is to identify which of them is the FJC model and which is the FRC model, from what you’ve learned in task 1.

Task 3) First-passage times of finding target

You shall now use the configuration x.npy (where x is either a or b)from task 2 that corresponds to the FRC, to set up another simulation, in which one bead of the polymer is of target type. Freely diffusing A particles have to find the target particle. The application you should have in mind is proteins docking to a certain part of nucleid acids, which is crucial for the function of each biological cell. To determine when an A particle has found the target we implement the following kind of reaction

$$ \mathrm{A} + \mathrm{target} \to \mathrm{B} + \mathrm{target} $$

The time when the first B particle is created, then corresponds to the first-passage time of that reaction.

3a) Perform a simulation that initializes the polymer from x.npy where the 10th bead is of type target, and places 50 A particles normally distributed (with variance $\sigma=0.5$) around the origin.

Therefore use the following system configuration

system = readdy.ReactionDiffusionSystem(
    box_size=[16., 16., 16.],
    periodic_boundary_conditions=[False, False, False],
    unit_system=None)

system.add_topology_species("monomer", 0.1)
system.add_topology_species("target", 0.1)
system.add_species("A", 0.5)
system.add_species("B", 0.5)

system.topologies.add_type("polymer")

origin = np.array([-7.5, -7.5, -7.5])
extent = np.array([15., 15., 15.])

system.potentials.add_box("monomer", force_constant=50., origin=origin, extent=extent)
system.potentials.add_box("target", force_constant=50., origin=origin, extent=extent)
system.potentials.add_box("A", force_constant=50., origin=origin, extent=extent)
system.potentials.add_box("B", force_constant=50., origin=origin, extent=extent)

with the following topology potentials

system.topologies.configure_harmonic_bond(
    "monomer", "monomer", force_constant=50., length=bond_length)
system.topologies.configure_harmonic_bond(
    "monomer", "target", force_constant=50., length=bond_length)


system.topologies.configure_harmonic_angle(
    "monomer", "monomer", "monomer", force_constant=20.,
    equilibrium_angle=2.530727415391778)
system.topologies.configure_harmonic_angle(
    "monomer", "monomer", "target", force_constant=20.,
    equilibrium_angle=2.530727415391778)
system.topologies.configure_harmonic_angle(
    "monomer", "target", "monomer", force_constant=20.,
    equilibrium_angle=2.530727415391778)

Define a boolean flag interaction = True. If the bool interaction is True then there should be a weak interaction between A and monomer particles with a force_constant of 50, desired_distance=bond_length, a depth of 1.4, and a cutoff of 2.2*bond_length. What does such an interaction result in?

Additionally there will be repulsion potentials between monomers and between A particles.

system.potentials.add_harmonic_repulsion(
    "monomer", "monomer", force_constant=50., 
    interaction_distance=1.1*bond_length)
system.potentials.add_harmonic_repulsion(
    "A", "A", force_constant=50., interaction_distance=1.5*bond_length)

Finally the system needs the reaction

system.reactions.add("found: A +(0.48) target -> B + target", rate=10000.)

where we use a very high rate, such that the reaction will happen directly on contact, where 0.48 is the contact distance.

Observe the number of B particles with a stride of 1.

Load the polymer configuration and turn the 10th monomer into a target.

init_pos = np.load("x.npy")

types = len(init_pos) * ["monomer"]
types[10] = "target" # define the target to be the 10-th monomer

top = sim.add_topology("polymer", types, init_pos)
for i in range(1, len(init_pos)):
    top.get_graph().add_edge(i - 1, i)

Place 50 A particles normally distributed, with variance $\sigma=0.5$, around the origin.

Simulate the system with a timestep of 0.001 for 30000 steps. Have a look at the VMD output.

  • How do the A particles interact with the polymer?
  • Calculate the first passage time from the observed number of particles, i.e. the time when the first B was created.

3b) Combine the simulation procedure above into a function of the signature

def find_target(interaction=False):
    ...
    # Since we will run many simulations
    # you may want to supress the textual output by
    # setting the two options show_progress
    # and show_summary
    sim.show_progress = False
    sim.run(..., show_summary=False)
    ...
    return passage_time

Hint: One such simulation should not take much longer than 10 seconds.

Gather passage times in a list

ts_int = []

Repeat the simulation many times (50-100 should suffice) and append the result to the list.

from tqdm import tqdm_notebook as tqdm
n=50
for _ in tqdm(range(n)):
    ts_int.append(find_target(interaction=True))

As this might take a while, you will want to observe how long the whole process takes, which is done here using tqdm.

Do the same for the case of no interaction, i.e. set interaction=False and gather the results in another list ts_noint.

ProTip: To save computation time, run this second case in a copy of your notebook (i.e. at the same time) and save the resulting list ts_noint into a pickle file, which you can read in in the original notebook.

For both cases interaction=True and interaction=False, calculate the distribution of first passage times, i.e. plot a histogram of the lists you constructed using plt.hist(). Use bins=np.logspace(0,2,20) and density=True.

When assuming a memory less (Poisson) process, the only relevant parameter is the mean rate $\lambda=N/\sum_{i=1}^N\tau_i$. Plot the distribution of first-passage-times with mean rate $\lambda$, i.e. the Poisson probability density (not the cumulative) of 1 event occurring before time t. Compare against your measured distribution of waiting times?

Is the process of finding the target with or without interaction well suited to be modeled as a memory-less process?

Now additionally mark the mean first passage time for each case using plt.vlines(). Is the difference of the two cases well described by the mean first passage time?

Thursday

This session will deal with the self assembly of macromolecular structures due to reactions.

Task 1) linear filament (e.g. actin) assembly

In this task we will look at a linear polymer chain, but instead of placing all beads initially, we will let it self-assemble from monomers in solution. The polymer that we will build will have the following structure

head--core--(...)--core--tail

where (...) means that there can be many core particles, but the structure is always linear.

The simulation box is periodic with box_size=[20., 20., 20.].

We define three topology particle species and one normal particle species

system.add_species("substrate", 0.1)
system.add_topology_species("head", 0.1)
system.add_topology_species("core", 0.1)
system.add_topology_species("tail", 0.1)

We also define one topology type

system.topologies.add_type("filament")

There should be the following potentials for topology particles

system.topologies.configure_harmonic_bond(
    "head", "core", force_constant=100, length=1.)
system.topologies.configure_harmonic_bond(
    "core", "core", force_constant=100, length=1.)
system.topologies.configure_harmonic_bond(
    "core", "tail", force_constant=100, length=1.)

The polymer should be rather stiff, which is the case for Actin filaments in biological cells. We can compactly write this for all triplets of particle types:

triplets = [
    ("head", "core", "core"),
    ("core", "core", "core"),
    ("core", "core", "tail"),
    ("head", "core", "tail")
]
for (t1, t2, t3) in triplets:
    system.topologies.configure_harmonic_angle(
        t1, t2, t3, force_constant=50., 
        equilibrium_angle=np.pi)

We now introduce a topology reaction. They allow changes to the graph of a topology in form of a reaction. Here we will use the following definition.

system.topologies.add_spatial_reaction(
    "attach: filament(head) + (substrate) -> filament(core--head)",
    rate=5.0, radius=1.5
)

Using the documentation, familiarize yourself, what this means. Do not hesitate to ask, since topology reactions can become quite a tricky concept!

Next create a simulation object. We want to observe the following

simulation.record_trajectory(stride=100)
simulation.observe.topologies(stride=100)

Add one filament topology to the simulation

init_top_pos = np.array([
    [ 1. ,0. ,0.],
    [ 0. ,0. ,0.],
    [-1. ,0. ,0.]
])
top = simulation.add_topology(
  "filament", ["head", "core", "tail"], init_top_pos)
top.get_graph().add_edge(0, 1)
top.get_graph().add_edge(1, 2)

Additionally we need substrate particles, that can attach themselves to the filament

n_substrate = 300
origin = np.array([-10., -10., -10.])
extent = np.array([20., 20., 20.])
init_pos = np.random.uniform(size=(n_substrate,3)) * extent + origin
simulation.add_particles("substrate", positions=init_pos)

Then, run the simulation

if os.path.exists(simulation.output_file):
    os.remove(simulation.output_file)
dt = 5e-3
simulation.run(400000, dt)

One important observable will be the length of the filament as a function of time. It can be obtained from the trajectory as follows:

times, topology_records = traj.read_observable_topologies()
chain_length = [ len(tops[0].particles) for tops in topology_records ]

The last line is a list comprehension. tops is a list of topologies for a given time step. Hence, tops[0] is the first (and in this case, the only) topology in the system. tops[0].particles is a list of particles belonging to this topology. Thus, its length yields the length of the filament.

1a) Have a look at the VMD output. Describe what happens? Additionally plot the length of the filament as a function of time. Note that you shall now plot the simulation time and not the time step indices, i.e. do the following

times = np.array(times) * dt

where dt is the time step size.

1b) Using your data of the filament-length, fit a function of the form

$$ f(t)=a(1-e^{-bt})+3 $$

to your data. You should use scipy.optimize.curve_fit to do so

import scipy.optimize as so

def func(t, a, b):
    return a*(1. - np.exp(-b * t)) + 3.

popt, pcov = so.curve_fit(func, times, chain_length)

print("popt", popt)
print("pcov", pcov)

f = lambda t: func(t, popt[0], popt[1])

plt.plot(times, chain_length, label="data")
plt.plot(times, f(times), label=r"fit $f(t)=a(1-e^{-bt})+3$")
plt.xlabel("Time")
plt.ylabel("Filament length")
plt.legend()
plt.show()

Question: Given the result of the fitting parameters

  • How large is the equilibration rate?
  • What will be the length of the filament for $t\to\infty$?

1c) We now introduce a disassembly reaction for the tail particle. This is done by adding the following to your system configuration.

def rate_function(topology):
    """
    if the topology has at least (head, core, tail)
    the tail shall be removed with a fixed probability per time
    """
    vertices = topology.get_graph().get_vertices()
    if len(vertices) > 3:
        return 0.05
    else:
        return 0.

def reaction_function(topology):
    """
    find the tail and remove it,
    and make the adjacent core particle the new tail
    """
    recipe = readdy.StructuralReactionRecipe(topology)
    vertices = topology.get_graph().get_vertices()

    tail_idx = None
    adjacent_core_idx = None
    for v in vertices:
        if topology.particle_type_of_vertex(v) == "tail":
            adjacent_core_idx = v.neighbors()[0].get().particle_index
            tail_idx = v.particle_index

    recipe.separate_vertex(tail_idx)
    recipe.change_particle_type(tail_idx, "substrate")
    recipe.change_particle_type(adjacent_core_idx, "tail")

    return recipe


system.topologies.add_structural_reaction(
    "detach",
    "filament",
    reaction_function=reaction_function, 
    rate_function=rate_function)

Familiarize yourself with this kind of structural topology reaction

Repeat the same analysis as before, and also observe your VMD output.

  • How large is the equilibration rate?
  • What will be the length of the filament for $t\to\infty$?

Task 2) Assembly of virus capsids

This task will suggest a model for the assembly of monomers into hexamers, and in the bonus task: the assembly of these hexamers into even larger superstructures.

First we need a model for one monomer, which should look as follows

It is essentially one bigger core particle with two sites attached. This is sometimes called a patchy particle, i.e. a particle with reaction patches. In the language of ReaDDy, this group of particles is a topology, let’s give it the name CA

system.topologies.add_type("CA")
system.add_topology_species("core", 0.1)
system.add_topology_species("site", 0.1)

The angle between the triplet site--core--site, shall be fixed to 120°. The sites shall react with the sites of other monomers to form a dimer that looks like

For this dimer, the four particles shall be confined to a single 2D plane, which we do using topology potentials. Summarizing all potentials that you need to configure for the particles:

  • harmonic bond between core--site with force constant 100 and length 1
  • harmonic bond between core--core with force constant 100 and length 2
  • harmonic bond between site--site with force constant 100 and length 0.1
  • harmonic angle between site--core--site with force constant 200 and equilibrium angle 120 degrees
  • harmonic angle between site--core--core with force constant 200 and equilibrium angle 120 degrees
  • harmonic angle between core--core--core with force constant 200 and equilibrium angle 120 degrees
  • dihedral between core--core--core--core with force constant 200, mutliplicity of 1 and equilibrium angle of 0
  • dihedral between site--core--core--core with force constant 200, mutliplicity of 1 and equilibrium angle of 0
  • dihedral between site--core--core--site with force constant 200, mutliplicity of 1 and equilibrium angle of 0
  • (normal) harmonic repulsion between core and core with force constant 80 and interaction distance 2

The formation of the dimer will be done using a spatial topology reaction

of the form

system.topologies.add_spatial_reaction(
    "attach: CA(site)+CA(site)->CA(site--site) [self=true]", rate=10., radius=0.4)

After such a reaction the connectivity looks like (...)--core--site--site--core--(...). In a second step after the reaction we want to get rid of the two site particles in the middle. This will be done using a structural topology reaction of the following kind (make sure that you understand what happens in these code snippets. Do ask, if you have problems understanding!). The first ingredient is the rate function for the structural reaction, i.e. given a toplogy this function shall return a very high rate, if there is a site--site connection, and shall return 0 otherwise

def clean_sites_rate_function(topology):
    edges = topology.get_graph().get_edges()
    vertices = topology.get_graph().get_vertices()

    if len(vertices) > 3:
        for e in edges:
            v1_ref, v2_ref = e[0], e[1]
            v1 = v1_ref.get()
            v2 = v2_ref.get()
            v1_type = topology.particle_type_of_vertex(v1)
            v2_type = topology.particle_type_of_vertex(v2)
            if v1_type == "site" and v2_type == "site":
                return 1e12
    else:
        return 0.
    return 0.

The second ingredient is the reaction function, that performs the removing of the two site particles, after the rate function has returned a very high rate

def clean_sites_reaction_function(topology):

    recipe = readdy.StructuralReactionRecipe(topology)
    vertices = topology.get_graph().get_vertices()

    def search_configuration():
        # dfs for finding configuration core-site-site-core
        for v1 in vertices:
            if topology.particle_type_of_vertex(v1) == "core":
                for v2_ref in v1.neighbors():
                    v2 = v2_ref.get()
                    if topology.particle_type_of_vertex(v2) == "site":
                        for v3_ref in v2.neighbors():
                            v3 = v3_ref.get()
                            if v3.particle_index != v1.particle_index:
                                if topology.particle_type_of_vertex(v3) == "site":
                                    for v4_ref in v3.neighbors():
                                        v4 = v4_ref.get()
                                        if v4.particle_index != v2.particle_index:
                                            if topology.particle_type_of_vertex(v4) == "core":
                                                return v1.particle_index, v2.particle_index, v3.particle_index, v4.particle_index

    core1_p_idx, site1_p_idx, site2_p_idx, core2_p_idx = search_configuration()

    # find corresponding vertex indices from particle indices
    core1_v_idx = None
    site1_v_idx = None
    site2_v_idx = None
    core2_v_idx = None
    for i, v in enumerate(vertices):
        if v.particle_index == core1_p_idx and core1_v_idx is None:
            core1_v_idx = i
        elif v.particle_index == site1_p_idx and site1_v_idx is None:
            site1_v_idx = i
        elif v.particle_index == site2_p_idx and site2_v_idx is None:
            site2_v_idx = i
        elif v.particle_index == core2_p_idx and core2_v_idx is None:
            core2_v_idx = i
        else:
            pass

    if (core1_v_idx is not None) and (core2_v_idx is not None) and (site1_v_idx is not None) and (
            site2_v_idx is not None):
        recipe.add_edge(core1_v_idx, core2_v_idx)
        recipe.separate_vertex(site1_v_idx)
        recipe.separate_vertex(site2_v_idx)
        recipe.change_particle_type(site1_v_idx, "dummy")
        recipe.change_particle_type(site2_v_idx, "dummy")
    else:
        raise RuntimeError("core-site-site-core wasn't found")

    return recipe

Finally add the structural reaction to the system

system.topologies.add_structural_reaction(
    "clean_sites", topology_type="CA", 
    reaction_function=clean_sites_reaction_function,
    rate_function=clean_sites_rate_function,
    raise_if_invalid=True,
    expect_connected=False)

2a) Simulate the system described above in a periodic box of size [25, 25, 25] for n_steps=50000 steps with a timestep of 0.005. Initially place 150 CA patchy particles uniformly distributed in the box. While simulating, observe the trajectory and topologies with the same stride.

sim.record_trajectory(n_steps//2000)
sim.observe.topologies(n_steps//2000)

What do you observe in the VMD output? Do particles assemble in the way you expected?

Hints:

  • You should place the particles such that the site particles are already at their prescribed 120 degree angle and a distance of 1 away from the core particle. Adding one such particle can be done in the following way

    core = np.array([0., 0., 0.])
    site1 = np.array([0., 0., 1.])
    site2 = np.array([np.sin(np.pi * 60. / 180.), 0., - 1. * np.cos(np.pi * 60. / 180.)])
      
    top = sim.add_topology("CA", ["site", "core", "site"], np.array([site1, core, site2]))
    top.get_graph().add_edge(0, 1)
    top.get_graph().add_edge(1, 2)
    
  • To distribute the particles uniformly you should add a random translation vector to all positions core, site1 and site2.
  • If you want to be super cool, you can rotate the patchy particle by a random amount before translating it. Ask google how to generate a random rotation matrix and how to apply it to your vectors core, site1 and site2

2b) From your output file and using the topologies and trajectory observable, calculate the time dependent distribution of molecular mass. This means: Given an instance of a topology, the degree of polymerization is the number of connected core particles in this topology. For one polymer the molecular mass is equal to the degree of polymerization. Obtain such a value for all topologies in a given timestep and make a histogram of that. Now that histogram only counts the occurrence of how many times a topology with a certain molecular mass shows up. To convert that into a distribution of molecular mass itself, you have to multiply the number of occurrence for each degree of polymerization by the degree of polymerization itself. Repeat this for all observed times to obtain as many histograms as there are timesteps.

Hints

  • The actual trajectory can be obtained from the trajectory file like so

    traj_file = readdy.Trajectory(out_file)
    traj = traj_file.read()
    
  • The particle type (string) of a particle with index v at time t is

    traj[t][v].type
    
  • Construct the histogram for each time using np.histogram(current_sizes, bins=bin_edges), where current_sizes is the list of the molecular masses you have obtained, and bin_edges=np.arange(1,10,1)
  • For plotting it might come in handy to convert the bin_edges to bin_centers, by calculating the midpoints for each bin
  • Plot the histograms using the following snippet

    xs = np.array(times) * dt
    ys = bin_centers
    X, Y = np.meshgrid(xs, ys-1)
    Z = all_histograms.transpose()
    plt.pcolor(X, Y, Z, cmap=plt.cm.viridis_r)
    plt.xlabel("Time")
    plt.ylabel("Degree of polymerization")
    plt.title("Distribution of molecular mass")
    

2c) Calculate a similar distribution of molecular mass, but now only for completely assembled topologies, i.e. topologies with no open site particles left. What is the percentage of “misfolded” topologies?

2d) From looking at your distribution of pentamers, hexamers, and heptamers. Can you form a full capsid out of the patchy particles we have used? Have a look at the introductory image with the viruses, what do you notice about the capsomers?

2e) Bonus task: Introduce a third reactive patch for each patchy particle called offsite, which allows binding to other offsite particles. In this way try to assemble a larger super structure out of the hexamers.

Friday

The idea for the friday session is to either

  • Continue to create spatial patterns with the diffusion-limited LV model
  • Try to simulate a system that you are working on / you know well
  • Think of any other scientifically related/unrelated system

At ca. 15:30h everyone will present their system: a short vmd movie and an explanation of what your system does, how it works.

Cookbook

This section has some solutions to common problems you might run into.

  • The order in which you create and manipulate the system and simulation matters!

Remember that the workflow should always be

system = readdy.ReactionDiffusionSystem(...)
# ... system configuration
simulation = system.simulation(...)
# ... simulation setup
simulation.run(...)
# ... analysis of results

If you made a mistake while registering species, potentials, etc., just create a new system and run the configuration again. The same goes for the simulation, just create a new one. In the jupyter notebooks it is sufficient, to just run all the cells again.

  • Simulation box and box potential for confining particles

If you want to confine particles to some cube without periodic boundaries, the box potential must be fully inside the simulation box. Remember that the box is centered around (0,0,0). For example the following will confine A particles to a cube of edge length 10

system = readdy.ReactionDiffusionSystem(
    box_size=[12., 12., 12.], unit_system=None,
    periodic_boundary_conditions=[False, False, False])
system.add_species("A", 1.)
origin = np.array([-5., -5., -5.])
extent = np.array([10., 10., 10.])
system.potentials.add_box(
    "A", force_constant=10.,
    origin=origin, extent=extent)

The two vectors origin and extent span a cube in 3D space, the picture you should have in mind is the following

  • Initial placement of particles inside a certain volume

If you have already defined origin and extent of your confining box potential, it is easy to generate random positions uniformly in this volume

n_particles = 30
uniform = np.random.uniform(size=(n_particles,3))
init_pos = uniform * extent + origin

Here uniform is a matrix Nx3 where each row is a vector in the unit cube $\in{[0,1]\times[0,1]\times[0,1]}$ . This is multiplied with extent, yielding a uniform cube ${[0,\mathrm{extent}_0]\times[0,\mathrm{extent}_1]\times[0,\mathrm{extent}_2]}$. If you add the origin to this you get this cube at the right position (with respect to our box coordinates centered around (0,0,0)), i.e.

$$ \begin{aligned} \{&[\mathrm{origin}_0,\mathrm{extent}_0+\mathrm{origin}_0]\\ \times&[\mathrm{origin}_1,\mathrm{extent}_1+\mathrm{origin}_1]\\ \times&[\mathrm{origin}_2,\mathrm{extent}_2+\mathrm{origin}_2]\} \end{aligned} $$
  • 2D plane

If you want to confine particles to a 2D plane, just use the box potential but make the extent in one dimension very small, i.e.

system = readdy.ReactionDiffusionSystem(
    box_size=[12., 12., 3.], unit_system=None,
    periodic_boundary_conditions=[False, False, False])
system.add_species("A", 1.)
origin = np.array([-5., -5., -0.01])
extent = np.array([10., 10., 0.02])
system.potentials.add_box(
    "A", force_constant=10.,
    origin=origin, extent=extent)

Having defined origin and extent it is now easy to add particles to this 2D plane. Note that I also made the box_size in z direction smaller, however it should be large enough.

  • Output file size

Please make use of your /storage/mi/user directories, your home will fill up quicker.

data_dir = "/storage/mi/user" # replace user with your username
simulation.output_file = os.path.join(data_dir, "myfile.h5")

Additionally, use a stride on your observables, e.g.

simulation.record_trajectory(stride=100)
# or
simulation.observe.particle_positions(stride=100)
  • Reaction descriptor language

In expressions like

system.reactions.add("fus: A +(2) B-> C", rate=0.1)

the value in the parentheses +(2) is the reaction distance.

  • look at VMD output if something is fishy

This might reveal some obvious mistakes. Therefore you must have registered the according observable

simulation.record_trajectory(stride=100) # use appropriate stride
# ... run simulation
traj = readdy.Trajectory(simulation.output_file)
traj.convert_to_xyz(particle_radii={"A": 0.5, "B": 1.})
# particle_radii is optional here

Then in a bash shell do

vmd -e myfile.h5.xyz.tcl

or prefix with ! in the jupyter notebook.

© Copyright 2020 AI4Science (former CMB) Group