Agent-based Spatiotemporal Molecular Distributions Evolving Under Simulation

Agent-based Spatiotemporal MOlecular Distributions Evolving Under Simulation, or ASMODEUS, is an R package for agent-based modelling of protein aggregation dynamics on locally flat surfaces, with an emphasis on simulating organisation of transmembrane receptors. We use a novel supervised learning approach, which inherits spatial statistics from existing data sets, such as those derived from single molecule localisation microscopy (SMLM).

Installation

The package can be installed using pre-existing functions in the devtools package. It is highly recommended that you use RStudio. If you do not have devtools installed, enter the following line into the command prompt: {r} install.packages("devtools") Then, simply enter the line below into the command prompt: {r} devtools::install_github("lucapanconi/asmodeus")

Usage

Seting Up

All code featured on this page is summarised in the test_suite.R file.

First, load the package into your R session. ```{r}

Implement ASMODEUS library.

library(asmodeus) ```

You may also set a seed for reproducibility.

```{r}

Set seed.

set.seed(1) ```

It is often convenient to set a main file directory for saving results. Enter a file directory of your choice here.

```{r}

Define main save location.

main <- "C:/Users/.../" ```

Static Data

ASMODEUS offers built-in functions for simulating static point clouds with circular clusters. These may be used to provide target statistics and test workflows.

```{r}

Simulate a data set with small circular clusters.

smallclusterdata <- simulatecircularclusters(ROI = 1000, numberofclusters = 10, clusterradius = 20, pointspercluster = 15, backgroundtoclusterratio = 0.33)

Simulate a data set with large circular clusters.

bigclusterdata <- simulatecircularclusters(ROI = 1000, numberofclusters = 5, clusterradius = 40, pointspercluster = 30, backgroundtoclusterratio = 0.33)

Simulate a data set with very large circular clusters.

verybigclusterdata <- simulatecircularclusters(ROI = 1000, numberofclusters = 3, clusterradius = 80, pointspercluster = 50, backgroundtocluster_ratio = 0.33) ```

We may also simulate complete spatially random data sets.

```{r}

Simulate a data set with complete spatially random (CSR) data. This will have the same number of points as the cluster data.

csrdata <- simulatecsr(ROI = 1000, numberofpoints = nrow(smallcluster_data)) ```

We can plot static point clouds with the plotpoints function.

```{r}

We can plot these point clouds individually.

plotpoints(smallclusterdata, xlimits = c(0, 1000), ylimits = c(0, 1000)) plotpoints(bigclusterdata, xlimits = c(0, 1000), ylimits = c(0, 1000)) plotpoints(verybigclusterdata, xlimits = c(0, 1000), ylimits = c(0, 1000)) plotpoints(csrdata, xlimits = c(0, 1000), ylimits = c(0, 1000)) ```

We can extract the Ripley's K function from any static point cloud with the ripleyk function.

```{r}

Get Ripley's K functions.

smallclusterK <- ripleyk(points = smallclusterdata, ROI = 1000) bigclusterK <- ripleyk(points = bigclusterdata, ROI = 1000) verybigclusterK <- ripleyk(points = verybigclusterdata, ROI = 1000) csrk <- ripleyk(points = csrdata, ROI = 1000) ```

These can be plotted separately or altogether using the plotlines function.

```{r}

Plot Ripley's K functions.

plotlines(data = cbind(smallclusterK, bigclusterK[,2], verybigclusterK[,2], csrk[,2]), labels = c("r", "K(r)")) ```

Basic Dynamic Simulations

The primary function of ASMODEUS is to simulate agent-based models of protein aggregation dynamics in which each agent (molecule) is iteratively fit to a target spatial statistic (the Ripley's K function). The basic, single population simulation initialises as a CSR distribution and fits to the global K function of a given target point cloud.

```{r}

Create a simulation with default parameters which tries to fit an initially CSR point cloud to the small cluster data set.

sim1 <- globalripleysim(pointcloud = smallcluster_data) ```

Here, we try to rederive the smallclusterdata point clouds from Markov Chain Monte Carlo simulation. We can animate this simulation using the animate_points function.

```{r}

Animate point distribution.

animate_points(simulation = sim1)

Save the animation using gganimate's anim_save function.

gganimate::animsave(filename = paste(main, "1SmallClustersSimulation.gif", sep = "")) ```

If we know the target K function, we can animate the change in the global Ripley's K function across the simulation.

```{r}

Animate change in Ripley's K function.

animatetarget(simulation = sim1, ROI = 1000, target = smallcluster_K)

We can remove the legend by setting legendposition = "none". The target here is always static anyway.

animatetarget(simulation = sim1, ROI = 1000, target = smallclusterK, legendposition = "none") gganimate::animsave(filename = paste(main, "2SmallClusters_K.gif", sep = "")) ```

Alternatively, we can animate both the K and H function by directly plugging in the target point cloud. The K function appears as above.

```{r}

Animate K function directly from point cloud (without manually calculating K function).

animateKfrompoints(simulation = sim1, ROI = 1000, pointclouds = smallclusterdata, legendposition = "none")

Animate H function directly from point cloud (without manually calculating H function).

animateHfrompoints(simulation = sim1, ROI = 1000, pointclouds = smallclusterdata, legendposition = "none") gganimate::animsave(filename = paste(main, "3SmallClustersH.gif", sep = "")) ```

We can define simulation parameters manually, use help("globalripleysim") for details. In particular, we can specify an initial distribution to start the simulation from. In this case, we have set the initial distribution as the big clusters data.

```{r}

Define simulation parameters manually, use help("globalripleysim") for details. In particular, we can specify an initial distribution to start the simulation from. In this case, we have set this initial distribution as the big clusters data.

sim2 <- globalripleysim(Dmin = 0, Dmax = 63, ROI = 1000, times = 1200, pointcloud = smallclusterdata, rmax = 200, nrval = 20, initialdistribution = bigclusterdata)

The simulation now starts from large clusters.

animatepoints(simulation = sim2) gganimate::animsave(filename = paste(main, "4BigClusterstoSmallClustersSimulation.gif", sep = "")) ```

This change of initial distribution is reflected in the H function.

```{r}

This is reflected in the H function.

animateHfrompoints(simulation = sim2, ROI = 1000, pointclouds = smallclusterdata, legendposition = "none") gganimate::animsave(filename = paste(main, "5BigClusterstoSmallClusters_H.gif", sep = "")) ```

If we do not have a target point cloud, but do have a target statistic, we can input this directly. Remember to define the number of points to be used in the simulation.

```{r}

In some cases, we can define our own target K function too, but we must specify how many points should be used in the simulation.

sim3 <- globalripleysim(Dmin = 0, Dmax = 63, ROI = 1000, times = 1200, target = smallclusterK, numberofpoints = 982, rmax = 200, nrval = 20, initialdistribution = bigcluster_data)

In this case, we still try to fit to the small clusters K function, so there is no difference.

animatepoints(simulation = sim3) gganimate::animsave(filename = paste(main, "6BigClusterstoSmallClustersSimulationfromTarget.gif", sep = "")) ```

The target used here is defined manually.

```{r}

The target here is defined manually.

animateHfrompoints(simulation = sim3, ROI = 1000, pointclouds = smallclusterdata, legendposition = "none") gganimate::animsave(filename = paste(main, "7BigClusterstoSmallClustersHfrom_Target.gif", sep = "")) ```

Multiple Targets

We can provide multiple possible target statistics or point clouds. These may be considered simultaneously or at distinct time points. For distinct time points, we specify the times for which the target must change.

```{r}

We can define multiple targets by inputting a list of point clouds or targets.

targetclouds <- list(bigclusterdata, smallcluster_data)

We can force the simulator to switch from one target to another at given time points. This can be done for as many targets as we like.

sim4 <- distripleysim(pointclouds = targetclouds, time_points = c(500, 1000))

The target is big clusters up until time point 500, then switches to small clusters until the end of the simulation (time point 1000).

animatepoints(simulation = sim4) gganimate::animsave(filename = paste(main, "8DiscreteTime_Points.gif", sep = "")) ```

When animating the global H function, we can observe the shift in the target.

```{r}

We can observe the shift in target in the H function. Make sure to specify the time points in the animation function.

animateHfrompoints(simulation = sim4, ROI = 1000, pointclouds = targetclouds, timepoints = c(500, 1000), legendposition = "none") gganimate::animsave(filename = paste(main, "9DiscreteTimePoints_H.gif", sep = "")) ```

When providing targets simultaneously, each agent will align with whichever target is closest. We can introduce perturbations at given times which offset each point randomly by a given strength (distance).

```{r}

We can also provide both targets simultaneously. Each point will move towards which ever target is closest.

sim5 <- simripleysim(pointclouds = targetclouds, times = 1000, perturbationtimes = 500, perturbationstrength = 100)

Points may drift towards either target. We can introduce perturbations which offest each point with distance given by perturbation_strength. Here we use a 100nm offset at time point 500.

animatepoints(simulation = sim5) gganimate::animsave(filename = paste(main, "10SimultaneousTargets.gif", sep = "")) ```

The impact of simultaneous targets is often unpredictable.

```{r}

The H function may assume one target over another, or fall somewhere inbetween.

animateHfrompoints(simulation = sim5, ROI = 1000, pointclouds = targetclouds, legendposition = "none") gganimate::animsave(filename = paste(main, "11SimultaneousTargets_H.gif", sep = "")) ```

We may also perturb points by cross-linking. This gives the effect of "squeezing" the points into clusters.

```{r}

Instead of offsetting points, we can cross-link points together. This gives the effect of "squeezing" points into clusters.

sim6 <- squeezeripleysim(pointclouds = targetclouds, times = 1000, perturbationtime = 500, numbercrosslinkedclusters = 5, numbercrosslinked = 30, crosslinkedcluster_radius = 40)

We form 5 cross-linked clusters at time point 500. Each cluster contains 30 points and is of radius 40.

animatepoints(simulation = sim6) gganimate::animsave(filename = paste(main, "12SimultaneousTargetswithSqueeze.gif", sep = ""))

A shift in the H function is shown at time point 500, when the squeeze perturbation is introduced.

animateHfrompoints(simulation = sim6, ROI = 1000, pointclouds = targetclouds, legendposition = "none") gganimate::animsave(filename = paste(main, "13SimultaneousTargetswithSqueeze_H.gif", sep = "")) ```

If we don't wish to immobilise agents after cross-linking, we can release them. This allows them to continue drifting.

```{r}

If we don't wish to immobilise points after cross-linking, we can release them.

sim7 <- squeezereleaseripleysim(pointclouds = targetclouds, times = 1000, perturbationtime = 500, numbercrosslinkedclusters = 5, numbercrosslinked = 30, crosslinkedclusterradius = 40) animatepoints(simulation = sim7) gganimate::animsave(filename = paste(main, "14SimultaneousTargetswithSqueezeandRelease.gif", sep = "")) animateHfrompoints(simulation = sim7, ROI = 1000, pointclouds = targetclouds, legendposition = "none") gganimate::animsave(filename = paste(main, "15SimultaneousTargetswithSqueezeandReleaseH.gif", sep = "")) ```

Alternatively, we can manually alter agent step sizes (speed) at distinct time points.

```{r}

We can manually change the step size (speed) of points during the simulation by inputting a vector of Dmin and Dmax values. We can also provide a perturbation to offset points.

sim8 <- speedripleysim(Dmin = c(0, 10), Dmax = c(10, 100), pointclouds = targetclouds, times = 1000, perturbationtimes = 500, perturbationstrength = 100) animatepoints(simulation = sim7) gganimate::animsave(filename = paste(main, "16SimultaneousTargetswithSpeedChange.gif", sep = "")) animateHfrompoints(simulation = sim7, ROI = 1000, pointclouds = targetclouds, legendposition = "none") gganimate::animsave(filename = paste(main, "17SimultaneousTargetswithSpeedChange_H.gif", sep = "")) ```

Multiple Populations

We can simulate multiple isolated or interacting molecular species by giving a target distribution for each. Populations simulated in isolation will not take account of other molecular species when calculating spatial statistics.

```{r}

We can simulate multiple populations representing different molecular species. Each population can be given its own target function.

targetclouds <- list(smallclusterdata, bigclusterdata, verybigclusterdata)

Here we simulate them independently of each other.

sim9 <- popripleysimisol(pointclouds = targetclouds)

When animating, we can specify labels for each population.

animatepopulations(simulation = sim9, labels = c("Small Clusters", "Big Clusters", "Very Big Clusters")) gganimate::animsave(filename = paste(main, "18MultiplePopulations_Isolated.gif", sep = "")) ```

Labels can be removed by removing the legend.

```{r}

We can also plot this without labels.

animatepopulations(simulation = sim9, legendposition = "none") gganimate::animsave(filename = paste(main, "19MultiplePopulationsIsolatedwithout_Legend.gif", sep = "")) ```

We can also simulate interactions between populations. In this case, each agent takes account of all other agents, regardless of their population type. To do this, we must convert the point clouds to a single data frame with three columns: x coordinate, y coordinate and population number.

```{r}

We can also allow different populations to interact. To do this, we must create a point cloud with three columns: x coordinate, y coordinate and population number.

allpoints <- rbind(cbind(smallclusterdata, 1), cbind(bigclusterdata, 2), cbind(verybigclusterdata, 3)) sim10 <- popripleysim(pointcloud = allpoints) animatepopulations(simulation = sim10, legendposition = "none") gganimate::animsave(filename = paste(main, "20MultiplePopulations.gif", sep = "")) ```

After this, we can induce and track interactions between populations. Here we assume the first point cloud is made of activators, the second is made of inhibitors, and the third is made of agents, which may be acted on by either activators or inhibitors. When an inactive agent (population 3) comes within 5nm of an activator (population 1), it becomes an active agent. Conversely, when an active agent (population 4) comes within 5nm of an inhibitor (population 2), it becomes an inactive agent.

```{r}

We can induce molecular interactions after the simulation is complete. Here we assume the first point cloud is made of activators, the second is made of inhibitors, and the third is made of agents, which may be acted on by either activators or inhibitors.

sim11 <- population_changer(simulation = sim10, A = c(3, 4), B = c(4, 3), C = c(1, 2), D = c(5, 5))

When an inactive agent (population 3) comes within 5nm of an activator (population 1), it becomes an active agent. Conversely, when an active agent (population 4) comes within 5nm of an inhibitor (population 2), it becomes an inactive agent.

animatepopulations(simulation = sim11, labels = c("Activator", "Inhibitor", "Inactive Agent", "Active Agent")) gganimate::animsave(filename = paste(main, "21MultiplePopulationswithInteractions.gif", sep = "")) ```

We can track the change in number of a particular population over time.

```{r}

We can track changes in a population over time.

changes <- track_changes(simulation = sim11, A = 4) plotlines(changes, labels = c("time", "number active")) ```

License

Licensed under GNU General Public License v3.0.