This tutorial illustrates how to solve numerically the approximate master equation for any binary state model running in any network with degree distribution
using DifferentialEquations
using LinearAlgebra
# Plotting libaries
using PyCall
using PyPlot
# Library with useful functions
include("./AME_functions.jl");
We assume a network topology generated with the configuration model given a degree distribution
function dist(z,k)
return exp(-z)*((big(z)^k)/factorial(big(k)))
end;
Since we need our differential equations to be finite, we set a range [$k_{\rm min}$,$k_{\rm max}$] and [$\tau_{\rm min}$,$\tau_{\rm max}$] of degrees and ages to consider in our numerical integration.
# Network description
z = 8
k_max = 20
k_min = 0
# Ages considered
tau_max = 200
tau_min = 0;
The initial state is set by the initial seed
The vector u0[1] incorporates the states
rho_0 = 10^-3
u0 = zeros(2,k_max+1,k_max+1,tau_max+1)
for k in k_min:k_max
for m in 0:k
# Susceptible
u0[1,k+1,m+1,1] = (1 - rho_0)*binomial_dist(k,m,rho_0)
# Infected
u0[2,k+1,m+1,1] = (rho_0)*binomial_dist(k,m,rho_0)
end
end
First, we define an activation probability that accounts for the aging mechanism
function pA(j)
return 1/(j+2)
end;
With
-
Infection probability:
$F(m/k - T) = p_A(j) \cdot \theta(m/k - T)$ . -
Reset probability:
$F_R (m/k - T) = p_A(j) \cdot (1 - \theta(m/k - T))$ . -
Aging probability:
$F_A (m/k - T) = 1 - p_A(j)$ .
# Threshold model with exogenous aging
# Model parameters:
T = 0.1
# Infection/recovery probabilities
F = zeros(k_max+1,k_max+1,tau_max+1)
R = zeros(k_max+1,k_max+1,tau_max+1)
# Reset probabilities
FR = zeros(k_max+1,k_max+1,tau_max+1)
RR = zeros(k_max+1,k_max+1,tau_max+1)
# Aging probabilities
FA = zeros(k_max+1,k_max+1,tau_max+1)
RA = zeros(k_max+1,k_max+1,tau_max+1)
for k in k_min:k_max
for m in 0:k
if m/(k+10^(-10)) >= T
for j in tau_min:tau_max
F[k+1,m+1,j+1] = pA(j) # To change state must activate and exceed threshold
end
end
if m/(k+10^(-10)) < T
for j in tau_min:tau_max
FR[k+1,m+1,j+1] = pA(j) # To reset must activate and do not exceed threshold
end
end
for j in tau_min:tau_max
FA[k+1,m+1,j+1] = 1 - pA(j) # To age must not activate at all
end
end
end
PARAMS = (F,R,FR,RR,FA,RA);
We integrate numerically using the tools of the library DifferentialEquations.jl. We define the time interval to integrate numerically tint, the integration time step dt, the saving time step ds and the integration solver.
tint = (0.0, 100.0)
dt = 0.01
ds = 0.1
# The integration method is a Runge-Kutta 4 method
method = RK4()
sol = solve_AME(tint,dt,ds,method,PARAMS);
We are interested on how the fraction of susceptible
I = infected_from(sol)
S = susceptible_from(sol)
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(15,5))
ax1.loglog(sol.t,I,label="I(t)",color="blue",linewidth=4,linestyle="-")
ax1.loglog(sol.t,S,label="S(t)",color="red",linewidth=4,linestyle="-")
ax1.loglog(sol.t,I .+ S,label="S(t) + I(t) = 1",color="green",linewidth=4,linestyle="--")
ax1.tick_params(labelsize=15)
ax1.legend(fontsize=15)
ax2.plot(sol.t,I,label="I(t)",color="blue",linewidth=4,linestyle="-")
ax2.plot(sol.t,S,label="S(t)",color="red",linewidth=4,linestyle="-")
ax2.plot(sol.t,I .+ S,label="S(t) + I(t) = 1",color="green",linewidth=4,linestyle="--")
ax2.tick_params(labelsize=15)
ax2.legend(fontsize=15)
PyObject <matplotlib.legend.Legend object at 0x7f9086127750>
Another interesting magnitudes that we can obtain from the numerical solution is the average interface density (fraction of S-I links in the network) and the persistence (fraction of nodes in the intial state)
int = interface_from(sol) # Interface calculation
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(15,5))
ax1.loglog(sol.t,int,linewidth=4,color="grey")
ax1.tick_params(labelsize=15)
ax1.set_xlabel("Time steps t",fontsize=15)
ax1.set_ylabel("Interface density",fontsize=15)
time, per = persistence_from(sol) # Persistence calculation
ax2.plot(time,per,linewidth=4,color="purple")
ax2.set_yscale("log")
ax2.tick_params(labelsize=15)
ax2.set_xlabel("Time steps t",fontsize=15)
ax2.set_ylabel("Persistence",fontsize=15);
-
To avoid problems, one shoud always choose a
$\tau_{\rm max}$ higher than the maximum time of evolution$t_{\rm max}$ to consider all possible states that agents can exhibit. -
The library DifferentialEquations.jl includes a large number of numerical solvers of differential equations. We used Runge-Kutta solvers but faster solvers may be found for different specific models.