Reaction-Diffusion Systems -------------------------- In reaction-diffusion systems, reaction processes are defined per location for a finite set of `N` locations that are connected via links in a network. Through these links, reactions between locations can take place. This means that if we assume a reactive system of `C` compartments per location, the system can be described using an instance of :class:`epipack.numeric_matrix_epi_models.MatrixEpiModel` with :math:`N\times C` compartments: Each compartment is identified by a location and its epidemiological description. .. note:: Usually, reaction diffusion equations are somewhat computationally heavy, which is why for research purposes you should consider other solutions or writing custom code. Nevertheless, ``epipack``'s philosophy of fast prototyping comes in handy, so if you just aim at finding the right model and visualize it in a small system, the following tutorial might help. Networks ======== One such example might be an SIR process taking place in different locations on earth which are connected via the air traffic network. In this case, the reaction equations read .. math:: S_u + I_u \stackrel{\eta}{\longrightarrow} 2I_u I_u \stackrel{\rho}{\longrightarrow} R_u X_u \stackrel{w_{uv}}{\longrightarrow} X_v Here, :math:`X_u` represents a compartment :math:`X\in\{S,I,R\}` at location :math:`u` with population size :math:`N_u`, and :math:`w_{uv}` represents a transition rate of an individual from location `u` to location `v`. In this picture, the total number of individuals traveling in the network remains constant and it is assumed that the population size per location is constant for all nodes, as well, if :math:`w_{uv}=w_{vu}`. The last equation describes an edge-centric random walk, which is associated with a uniform equilibrium concentration. Let's assume we have an undirected network :math:`A_{uv}` and the transition rate per link is given as :math:`w_{uv}=\gamma A_{uv}`. We want to define a chain here: .. code:: python from epipack import MatrixEpiModel N = 4 nodes = list(range(N)) links = [ (u,u+1,1.0) for u in range(N-1) ] base_compartments = list("SIR") compartments = [ (node, comp) for node in nodes for comp in base_compartments ] model = MatrixEpiModel(compartments) See that instead of using strings for compartments, we're using a tuple that contains the node and the compartment type. This works because tuples are `hashable`. Now we want to define our reaction processes .. code:: python infection_rate = 3 recovery_rate = 1 mobility_rate = 0.05 quadratic_processes = [] linear_processes = [] for node in nodes: quadratic_processes.append( ( (node, "S"), (node, "I"), infection_rate, (node, "I"), (node, "I") ), ) linear_processes.append( ( (node, "I"), recovery_rate, (node, "R") ) ) for u, v, w in links: for C in base_compartments: linear_processes.extend([ ( (u, C), w*mobility_rate, (v, C) ), ( (v, C), w*mobility_rate, (u, C) ), ]) See that due to symmetry, we have to set transition reactions for :math:`u\rightarrow v` as well as :math:`v\rightarrow u` for the same link. Finally, we pass the processes to the model to be converted .. code:: python model.set_processes(quadratic_processes+linear_processes) All that's left to do is to define initial conditions. Let's say we do this by assuming 20 percent of the population on the first node is infected. .. code:: python initial_conditions = { ( node, "S" ): 1.0 for node in nodes } initial_conditions[(nodes[0], "S")] = 0.8 initial_conditions[(nodes[0], "I")] = 0.2 model.set_initial_conditions(initial_conditions,allow_nonzero_column_sums=True) Note that we pass the ``allow_nonzero_column_sums=True`` flag to suppress a warning that the initial conditions do not sum up to unity. Now we can integrate and plot the ``I`` compartment for each node .. code:: python # set compartments for which you want to obtain the # result plot_compartments = [ (node, "I") for node in nodes ] # integrate import numpy as np t = np.linspace(0,12,1000) result = model.integrate(t,return_compartments=plot_compartments) # plot result import matplotlib.pyplot as plt plt.figure() for (node, _), concentration in result.items(): plt.plot(t, concentration, label=str(node)) plt.xlabel("time") plt.ylabel("I") plt.legend() plt.show() Finally: .. image:: vis_media/chain_I.png :width: 75% .. note:: As stated above, nodes are assumed to carry uniform equilibrium density of random walkers (as is the case in edge-centric random walks). If, instead, you want to assume that nodes have an equilibrium density proportional to their total in-/outflux (weighted degree), you have to rescale the links (which are then equal to transition probabilities and asymmetrical). For further info, see the supplementary material of https://science.sciencemag.org/content/342/6164/1337.figures-only Visualization ============= Reaction-diffusion systems are interesting to watch unfold. Because we know this, ``epipack`` provides a visualization function for reaction diffusion systems: :func:`epipack.vis.visualize_reaction_diffusion`. Let's create a modular hierarchical network from the cMHRN_ package and style it .. code:: python import networkx as nx import netwulf as nw import cMHRN # load edges from txt file and construct Graph object N, edges = cMHRN.fast_mhrn(8,3,7,0.18,True) G = nx.Graph() G.add_edges_from(edges) # visualize and save visualization network, config = nw.visualize(G) nw.save("MHRN.json",network,config) .. code:: python # load visualization network, config, _ = nw.load("MHRN.json") # get the network properties N = len(network['nodes']) nodes = list(range(N)) links = [ ( link['source'], link['target'], 1.0 ) for link in network['links'] ] Subsequently, we set up the model exactly as above. One thing that you should know is that internally, an instance of ``MatrixEpiModel`` creates a one-dimensional :math:`N\times C`-long vector that contains the state of each compartment. In order for the visualization function to know which entry it should plot for each node, we have to provide it with the appropriate compartments. To this end, we construct a list of `N` entries, each entry `i` maps node `i` to the compartment whose concentration is supposed to be shown for this node. Here, we want the node to be scaled according to its "I" compartment, which we expect to take values between ``0`` and ``0.3``: .. code:: python node_compartments = [ (node, "I") for node in nodes ] value_extent = [0,0.3] And finally, we can start the visualization .. code:: python from epipack.vis import visualize_reaction_diffusion dt = 0.04 visualize_reaction_diffusion(model, network, dt, node_compartments, value_extent=value_extent, ) .. video:: ../_static/reac_diff_MHRN.mp4 :width: 500 Lattices ======== Sometimes, a network is not the structure you want to be looking at. That's fine, you do you, go for the lattice. However, what is a lattice but a network? If you want so simulate/visualize on a lattice, simply do everything as described above, but use a lattice network .. code:: python from epipack.vis import visualize_reaction_diffusion, get_grid_layout from epipack import get_2D_lattice_network N_side = 30 N = N_side**2 nodes = range(N) links = get_2D_lattice_network(N) network = get_grid_layout(N) Now, set up everything as above and simulate on a lattice: .. code:: python visualize_reaction_diffusion(model, network, dt, node_compartments, value_extent=value_extent, config = { 'draw_nodes_as_rectangles': True, 'draw_links': False, } ) .. video:: ../_static/reac_diff_lattice.mp4 :width: 500 .. _cMHRN: https://github.com/benmaier/cMHRN