Edit on GitHub

Particle aggregation

In this model we are going to make a model with the following properties:

  • Soft spheres with repulsion and adeshion
  • Viscous brownian motion dynamics
  • Particle confinement in a simulation space

The model

The model will act with viscous brownian motion dynamics of the form:

\[dx_i = \underbrace{\sum_j f(x_i,x_j)}_{\text{deterministic term}} dt + \underbrace{D}_{\text{stochastic term}} dW\]

where the forces of interaction are of the form

\[f(x,y)=\begin{cases} F_{rep}(r_{rep}-x)\hspace{1cm}\text{if }||x-y||_2 <= r_{rep}\\ -F_{atr}(r_{atr}-x)\hspace{.75cm}\text{if }r_{rep} < ||x-y||_2 < r_{at}\\ 0\hspace{2.45cm}\text{otherwise }\\ \end{cases}\]

and some radius of repulsion $r_{rep}$ and attraction $r_{atr}$.

The forces of the particles over the walls will be of the maximum size of repulsion $f_{rep}$.

We are going to create the model step by step with the corresponding tests to check that it works properly.

Load the package and the plotting functions

using CellBasedModels
using GLMakie
Makie.inline!(true);

Construct the AgentBasedModel (ABM)

We are going to make the simulations of a model in 2D. By declaring a model with ABM(2,...) we are already declaring a model that includes two agent parameters: x and y. Agent parameters are parameters that have a value for each agent independently.

If we want to have track of the forces exerted over each particle, we need to define additionally two more agent parameters:

  • fx
  • fy

Additionally, our model has several global parameters. That is, parameters that are shared among all the agents in the system. Those are:

  • rRep: radius of repulsion
  • frep: force of repulsion
  • rAtr: radius of attraction
  • fAtr: force of atraction
  • D: diffusion coefficient

For defining the Stochastic Differential Equation, we will need to add two additional keyword arguments:

  • agentODE: Where we define the deterministic part of the SDE
  • agentSDE: Where we define the stochastic part of the SDE

Let's create the model.

model = ABM(2,

    #In the keyword argument agent we declare all the agent arguments and its scope
    agent = Dict(
            :fx => Float64,
            :fy => Float64,
        ),

    #In the keyword argument model we declare all the model arguments and its scope   
    model = Dict(
            :rRep => Float64,
            :fRep => Float64,
            :rAtr => Float64,
            :fAtr => Float64,
            :D => Float64
        ),

    # In here we define the ODE part
    agentODE = quote

        # Compute adhesion and repulsion forces
            #Reset the parameters to zero zero
        fx = 0
        fy = 0
            #Go over the neighbors and add the forces, for that we use the macro @loopOverNeighbors(iteratorSymbol,code) see more in documentation
        @loopOverNeighbors i2 begin
            d = CBMMetrics.euclidean(x,x[i2],y,y[i2])
            dirx = (x-x[i2])/d
            diry = (y-y[i2])/d
            if d < rRep #Repulsion forces
                fx += fRep*(rRep-d)*dirx  
                fy += fRep*(rRep-d)*diry  
            elseif d < rAtr #Attraction forces
                fx += -fAtr*(rAtr-d)*dirx  
                fy += -fAtr*(rAtr-d)*diry  
            end
        end
            #Add the forces comming from the boundary interaction
        if x < simBox[1,1]+rRep/2
            fx += fRep
        elseif x > simBox[1,2]-rRep/2
            fx -= fRep        
        end
        if y < simBox[2,1]+rRep/2
            fy += fRep
        elseif y > simBox[2,2]-rRep/2
            fy -= fRep
        end

        # Finally, define the deterministic term of the SDE
        dt(x) = fx
        dt(y) = fy
    end,

    # In here we define the SDE part
    agentSDE = quote
        # SDE term
        dt(x) = D
        dt(y) = D
    end
);

Tests

Now, we can test the correct definition of the properties. For that we are going to make two tests:

  • Test 1: Repulsion the agents
  • Test 2: Resulsion from boundaries

Test repulsion from boundaries

For making this test we are going to create a simulation with two agents, put them close together and see how they separate from each other.

For that we use the object Community(ABM, kwargs...) and we are going to define the following arguments:

  • N: The number of agents of the model.
  • dt: The integration step.
  • simBox: The simulation space.

Then, we will initialize the model and agent parameters.

simBox = [-5. 5; -2 2]
#Initialize community
com = Community(model,
        N=2,
        dt=.1,
        simBox = simBox
        );

#Setup user paramaters
        #Constants
com.rRep=.8
com.fRep=1
com.rAtr=1.
com.fAtr=1.
        #Agent parameters
com.x=[-.1,.1]
com.y=[0.,0.];
# If we do not define fx and fy, they are initialy set to zeros
2-element Vector{Float64}:
 0.0
 0.0

Now that everything is set, we can evolve the community with the function evolve! for some number of steps.

evolve!(com,steps=50)

And now, we visualize the results using Makie.

For that, we extract from the community all the parameters that we want to plot for all times.

d = getParameter(com,[:t,:x,:y,:fx,:fy])
Dict{Symbol, Vector} with 5 entries:
  :fy => [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0…
  :y  => [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0…
  :fx => [[-0.6, 0.6], [-0.6, 0.6], [-0.48, 0.48], [-0.36, 0.36], [-0.264, 0.26…
  :t  => [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0  …  4.1, 4.2, 4.3, 4…
  :x  => [[-0.1, 0.1], [-0.16, 0.16], [-0.22, 0.22], [-0.268, 0.268], [-0.304, …

And then we construct the plot.

#Make the plot
fig = Figure(resolution=(1100,550))

#Plot the spheres at different times
for (j,i) in enumerate(1:5:20) #Step times
    #Create axis
    ax = Axis(fig[1,j],title="t = $(round(d[:t][i],digits=2))")
    #Add a simulation box rectangle in the background
    mesh!(ax, [simBox[1,1] simBox[2,1];simBox[1,1] simBox[2,2];simBox[1,2] simBox[2,2];simBox[1,2] simBox[2,1]], [1 2 4; 2 3 4], color=:lightgrey)
    #Add spheres with attraction radius    
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rAtr][1]/2,color=(:blue,.3),transparency=true)
    #Add spheres with attraction radius
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rRep][1]/2,color=:red)
    #Define the limits
    xlims!(ax,-2,2)
    ylims!(ax,-2,2)
end

#Plot position over time of one of the particles
ax = Axis(fig[2,1:2],title="distance",xlabel="t",ylabel="|x₁-x₂|")
lines!(d[:t],[abs(j-i) for (i,j) in d[:x]])

#Plot forces over time of one of the particles
ax = Axis(fig[2,3:4],title="fx",xlabel="t",ylabel="|fx|")
lines!(d[:t],[abs(i) for (i,j) in d[:fx]])

#Show the plots
display(fig)

png

Test repulsion from boundaries

Now that we have shown the process, we do the same for the boundaries.

simBox = [-2. 2; -2 2]
#Initialize community
com = Community(model,
        N=1,
        dt=.1,
        simBox = simBox,
        );

#Setup user paramaters
com.rRep=.8
com.fRep=1
com.rAtr=1.
com.fAtr=1.
com.x=[0]
com.y=[-2.5];
#Evolve
evolve!(com,steps=15)
#Get parameters
d = getParameter(com,[:t,:x,:y,:fx,:fy])
Dict{Symbol, Vector} with 5 entries:
  :fy => [[1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0],…
  :y  => [[-2.5], [-2.4], [-2.3], [-2.2], [-2.1], [-2.0], [-1.9], [-1.8], [-1.7…
  :fx => [[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],…
  :t  => [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,…
  :x  => [[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],…
#Plot
fig = Figure(resolution=(1100,550))

#Plot visually the results
for (j,i) in enumerate(1:4:13)
    ax = Axis(fig[1,j],title="t = $(round(d[:t][i],digits=2))")
    mesh!(ax, [simBox[1,1] simBox[2,1];simBox[1,1] simBox[2,2];simBox[1,2] simBox[2,2];simBox[1,2] simBox[2,1]], [1 2 4; 2 3 4], color=:lightgrey)
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rAtr][1]/2,color=(:blue,.3),transparency=true)
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rRep][1]/2,color=:red)
    xlims!(ax,-3,3)
    ylims!(ax,-3,3)
end

#Plot position in i
ax = Axis(fig[2,1:2],title="distance",xlabel="t",ylabel="|x₁+2|")
lines!(d[:t],[i+2 for (i,) in d[:y]])

#Plot forces in i
ax = Axis(fig[2,3:4],title="fx",xlabel="t",ylabel="|fx|")
lines!(d[:t],[abs(i) for (i,) in d[:fy]])

display(fig)

png

Constructing over a predefined model

Now that we know that everything works for the model, we could start scaling it to many particles. But before we do that, we would like to add one more property to it as we would like to study the behavior of the particles as we cool down the movement.

For that, we want to make a new model that have exactly the same behaviour as before but in addition we want that the model parameter D that controld the diffusion strength is lowered down slowly.

We could rewrite all the model again from scratch, but in this case that is not necessary. We can define a new model that uses as a base the previous one and add the additional property.

Let's define the cooling down of the model with an exponential decay:

\[dD = -\alpha_TDdt\]

model2 = ABM(2,

    #Import the previous model
    baseModelInit = [model],

    #Add an additional model parameter that defines the speed of the diffusion parameter decay
    model = Dict(
        :αT => Float64
    ),

    #Add a ODE to describe the difussion
    modelODE = quote
        dt(D) = -αT*D
    end,

    agentAlg = CBMIntegrators.EM(),
    modelAlg = DifferentialEquations.Euler(),
    neighborsAlg = CBMNeighbors.CellLinked(cellEdge=2),
);

Evolve the model an check the behaviour

We can define now a now Commmunity and see the evolution in time. We can check that as the diffusion constant is reduced the particles start to aggregate as the adhesion forces dominate the behaviour over the active term.

simBox = [-10. 10; -10 10]
N = 100
com = Community(model2,
        N=N,
        dt=0.1,
        simBox = simBox,
        )

com.rRep=.8
com.fRep=1
com.rAtr=1.
com.fAtr=1.
com.αT=.001 
com.x = 20 .*rand(N) .-10
com.y = 20 .*rand(N) .-10
com.D=.15
0.15
evolve!(com,steps=8000,saveEach=1)
d = getParameter(com,[:x,:y])

fig = Figure(resolution=(1400,200))

for (pos,i) in enumerate(1:999:8000)
    ax = Axis(fig[1,pos],title="D=$(round(com[i].D[1],digits=4))")
    mesh!(ax, [simBox[1,1] simBox[2,1];simBox[1,1] simBox[2,2];simBox[1,2] simBox[2,2];simBox[1,2] simBox[2,1]], [1 2 4; 2 3 4], color=:lightgrey)
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rAtr][1]/2,color=(:blue,.0),transparency=true)
    meshscatter!(ax,d[:x][i],d[:y][i],markersize=com[:rRep][1]/2,
                color=[com[i].D[1] for i in 1:100],
                colorrange=(0.,0.15),
                colormap=:reds
                )
    xlims!(ax,(simBox[1,:].+[-1,1])...)
    ylims!(ax,(simBox[2,:].+[-1,1])...)
    # lines!(ax,[i[1] for i in d[:x]],[i[1] for i in d[:y]])
end

# save("agg.png",fig)

display(fig)

png

Neighbors computational dependence

It may be tempting to run simulations with a large number of particles. However, as we increase the number of particles, the computational time increases quadratically with the system size.

This is due to the fact that the interaction forces in @loopOverNeighbors check all the forces over all neighbors. However, in this model only particles that are close together really interact.

In CellBasedModels there are implementations of more efficient neighbors searching methods. You can provide them when creating the Community model.

Check how alternative algorithms improve the speeding time.

modelFull = ABM(2, baseModelInit=[model2], agentAlg = CBMIntegrators.EM(), neighborsAlg=CBMNeighbors.Full())
modelVerlet = ABM(2, baseModelInit=[model2], agentAlg = CBMIntegrators.EM(), neighborsAlg=CBMNeighbors.VerletDisplacement(skin=2,nMaxNeighbors=20))
modelCellLinked = ABM(2, baseModelInit=[model2], agentAlg = CBMIntegrators.EM(), neighborsAlg=CBMNeighbors.CellLinked(cellEdge=2))

function initialize(model,N,simBox)
        return Community(model,N=N,
                simBox = simBox,
                
                rRep=.8,
                fRep=1,
                rAtr=1.,
                fAtr=1.,
                D = .5,
                x=rand(N).*(simBox[1,2]-simBox[1,1]).+simBox[1,1],
                y=rand(N).*(simBox[2,2]-simBox[2,1]).+simBox[2,1],
                dt=.1
                );
end;
ρ=.1
S = 1

N = []
tFull = []
tVerlet = []
tCellLinked = []
for S in 1:.5:7
    simBox = S.*[-10. 10; -10 10]
    n = round(Int64,ρ*(simBox[1,2]-simBox[1,1])*(simBox[2,2]-simBox[2,1]))
    push!(N,n)

    com = initialize(modelFull,n,simBox); #Full (default) algorithm
    evolve!(com,steps=10,saveEach=10,saveCurrentState=true)
    t = @elapsed evolve!(com,steps=1000,saveEach=10,saveCurrentState=true)
    push!(tFull,t)

    com2 = initialize(modelVerlet,n,simBox); #Verlet list algorithm
    evolve!(com2,steps=10,saveEach=10,saveCurrentState=true)
    t = @elapsed evolve!(com2,steps=1000,saveEach=10,saveCurrentState=true)
    push!(tVerlet,t)

    com3 = initialize(modelCellLinked,n,simBox); #Cell linked algorithm
    evolve!(com3,steps=10,saveEach=10,saveCurrentState=true)
    t = @elapsed evolve!(com3,steps=1000,saveEach=10,saveCurrentState=true)
    push!(tCellLinked,t)
end
fig = Figure(resolution=(800,600))
ax = Axis(fig[1,1],xlabel="Number of agents",ylabel="seconds")

l1 = scatter!(ax,Float64.(N),Float64.(tFull))
l2 = scatter!(ax,Float64.(N),Float64.(tVerlet))
l3 = scatter!(ax,Float64.(N),Float64.(tCellLinked))
Legend(fig[1,2],[l1,l2,l3],["Full","Verlet List","Cell Linked"])

display(fig)

png