How host mobility formulations shape estimates of pathogen dispersal and epidemic risk in non endemic regions

The repository provides all scripts necessary to reproduce the numerical analyses and figures reported in the paper Presigny et al. (2026), which includes the:

• definition of the force of infection (FoI) and the diffusion model;
• generation of the synthetic networks;
• generation of the empirical italian networks at NUTS2 and NUTS3 level;
• computation of the observables (e.g. basic reproduction number, attack rate)
• estimation the mobility parameter epsilon from empirical mobility data;
• execution of the deterministic and Monte-Carlo epidemiological simulations;
• reproduction of all figures presented in the manuscript.

---

Installation

Clone the repository and install in editable mode:

git clone https://github.com/CoMuNeLab/Host_mobility_formulations
cd Host_mobility_formulations
pip install -e .

Package Structure

src/analysis/ \
├── epidemics.py       # Run deterministic epidemics \
├── mc_analysis.py   # Monte-Carlo simulations utility \
├── metrics.py           # Computation R0 and attack rate \

src/models/ \
├── CLASS_Allfield_model_cluster.py      # Run deterministic force-of-infection epidemics \
├── CLASS_Allfield_model_MC.py  # Run Monte-Carlo force-of-infection epidemics (fixed importation rate) \
├── CLASS_DIffusion_model_cluster.py           # Run deterministic diffusion epidemics \
├── CLASS_DIffusion_model_MC.py # Run Monte-Carlo diffusion epidemics (fixed importation rate) \

src/synthetic/ \
├── radiation_model.py       # Analytical radiation model \
├── soneira_peebles.py   # Build Soneira-Peebles node distribution \

src/utils/ \
├── io.py       # Input-output utilities \
├── maths.py   # Maths operation utilities \

Requirements

The code requires Python >=3.10.

Required packages include

• NumPy
• SciPy
• pandas
• geopandas
• matplotlib
• networkx

---

Data Availability

All open access raw data, and intermediate/processed data that the authors can publicly distribute are available in the Zenodo repository: https://doi.org/10.5281/zenodo.21355910

Intermediate results are stored to avoid recomputation whenever possible.

Description of the model classes

The models presented in the manuscript are encoded into classes where the dynamics is runned and stored. Refer to the manuscript and supplementary for the values of the parameters.

Input for the classes

• Y0: matrix (Nx7), initial state of the system. Each row represents the state of a node, each colum represents the numerosity of the different compartments:
  • 0:Sh;1:Eh;2:Ih;3:Rh;4:Sv;5:Ev;6:Iv
  • obtained in paraemter files by providing the population distristibution of the system, of the infection reservoir (Y0_W) and the vector-to-human ratio (M)
• parameter: dictionnary, contains all the parameters for the dynamical system to work.
  • "N_W_system": daily host arrival from the infection reservoir to the system (list)
  • "p0": average number of vectors transported per moving host
  • "eta": inverse of the characteristic length of the autonomous vector dispersal
  • "c_vh": effective contact rate vector to host (list)
  • "c_hv": effective contact rate host to vector (list)
  • "T_id": host infection duration (list)
  • "T_iit": intrinsic incubation time (list)
  • "T_lh": host life expectancy (list)
  • "K": vector carrying capacity (list)
  • "T_eit": extrinsic incubation time (list)
  • "T_lv": vector life expectancy (list)
  • "epsilon": time-scale separation parameter
  • "distance_matrix": distances between each patch in the system (in km) (matrix)
  • "mobility_matrix": conditional probability that an individual departing from one node to another node (matrix)
  • "T_W": proportion of hosts arrival from infection reservoir in node i
  • time_span: time of simulation
  • l_t: if l_t is specified, specify at which time-step the dynamics is saved
  • leap: if Montecarlo simulations. If None Gillespie algorithm, If float, give the leaping time duration (Tau-leaping algorithm)
  • N_sim = 12 if Montecarlo simulations indicate the number of simulations to run
• Model: Always fixed to the value = Mobility_Matrix_Model

Run simulations

1. Import the right model file (e.g. CLASS_DIffusion_model_cluster)
2. Initialize the instance: epidemics = Metapopulation(Y0, parameter, Model)
3. Run the mea-field model: epidemics.run_model(time_span,l_t)
4. Run the Monte-Carlo model: epidemics.run_montecarlo(N_sim,l_t,leap=leap)

Output for the classes

Once a simulationn has been performed with an instanciated class, the instance keep track of the dynamics of every comaprtment in every nodes. Here are the methods to get information from the class (we detail here for deterministic simulations, for Monte-Carlo simulations add the suffix "_Montecarlo" to the getter):

• __get_number_host_ODE(which_point): For the node "which_point", return the total number of hosts in function of time

• __get_number_vector_ODE(self,which_point): For the node "which_point", return the total number of vectors in function of time

• get_variables(self,which_var,which_point=0): For the node "which_point", return the toal number of host/vector in a given compartment in function of time (0:Sh;1:Eh;2:Ih;3:Rh;4:Sv;5:Ev;6:Iv)

• self.l_t: all the time step at which the dynamics was recorded

• self.Y_t: The dynamics of every node and every comaprtment in function of time

Reproducing the Results

The analyses are organized according to the figures presented in the manuscript. Please, refer to the Methods and the Supplementary Material of the article for more details.

Estimation of the mobility parameter

The estimation procedure reproduces the calibration of the mobility parameter epsilon from empirical mobility and population data. The workflow is described for the italian provinces but is equivalent for italian regions:

• create_mobility_matrix_IT_provinces.py :From colocation data, create the weekly mobility matrices.
• create_weekly_IT_population_provinces.py: From colocation data, create the population matrix for each week.
• find_epsilon_proovinces.py: Compute the value for epsilon.

Running this workflow produces the calibrated values used throughout the manuscript.

Construction of the italian systems (Italy_systems)

This directory contains the scripts to build the italian mobiltiy networks. The workflow is described for the italian provinces but is equivalent for italian regions:

• make_provinces_IT.py : Produce and save the position of the patches, the distance between them, their population, their name.
• mobility_provinces_IT.py : Produce the mobility matrix.
• airport_provinces_IT.py : Associate the nodes with airports and their relative weights with daily arrival of passengers.

Figure 2

This analysis reproduces the synthetic benchmark presented in Figure 2.

• make_synthetic_data : scripts to generate synthetic mobility networks, coarse-grained them and visualize them.
• compute_R0_and_attack_rate : Run the epidemics, compute attack rate and the basic reproduction number and parameter files for the synthetic networks.
  • parameters_R0_0.98.py : parameter file for the numerical experiments.
  • parameters_R0_0.98Nw=25000.py : parameter file for the numerical experiments with the inclusion of the infection reservoir.
  • run_epsilon_p0.py : Run deterministic simulations by varying p0 and epsilon
  • compute attack_rate_cluster.py_ : From the numerical experiments, compute attack rate
  • compute_numerical_R0_epsilon_cluster.py : Compute the basic reproduction number using both models while varying epsilon and p0.
  • assemble_R0.py : Assemble the results of the previous script into an unique file per networks.
• phase_diagram_display_epsilon.py : Show for both models the metrics (basic reproduction number and attack rate) in function of epsilon and p0

variables
var  # "d0=0.150R0=0.98" OR "d0=0.150R0=0.98Nw=25000", if inclusion of infection reservoir
which_model  # allfield (force-of-infection) OR diffusion, selection of model
type_network  # always fix to "random"
parameter_file # No infection reservoir: parameters_R0=0.98.py, with infection reservoir: parameters_R0=0.98Nw=25000.py
airport # True if display results with infection reservoir, False otherwise.
measure # R0, if display R0 results, "attack" if display attack rate results.
make_avg # if True, display the average of the selected measure over all the synthetic networks.
make_std à If True, display the standard deviation of the selected measure over all the synthetic networks.

Figures 3 and 4

These analyses reproduce the experiments performed on empirical mobility systems at different spatial resolutions. The workflow is described for the italian provinces but is equivalent for italian regions:

• compute_R0_and_attack_rate : Idem as Fig2
• phase_diagram_display_epsilon_Italy.py : Idem as Fig2

variables
var  # "d0=0.150R0=0.98" OR "d0=0.150R0=0.98_airport", if inclusion of infection reservoir
which_model  # allfield (force-of-infection) OR diffusion, selection of model
type_net  # regions OR provinces, selection of mobiltiy network
parameter_file # No infection reservoir: Fig5/compute_R0_and_attack_rate/parameters_R0=0.98_italy_noairport.py with infection reservoir: Fig5/compute_R0_and_attack_rate/parameters_R0=0.98_italy.py
airport # True if display results with infection reservoir, False otherwise.
measure # R0, if display R0 results, "attack" if display attack rate results.

Figure 5

This workflow reproduces the additional numerical experiments presented in Figure 5 using the same computational framework with the modified parameter settings described in the manuscript:

• compute_equilibrium_distribution.py : Produce the italian system where the population is fixed to the one of the diffusion model at equilibrium.
• compute_R0_and_attack_rate : Idem as Fig2

Variables are the same as in Figures 3 and 4.

Figures 6 and 7

This analysis reproduces the analysis presented in Figure 6 and Figure 7.

• analysis : the directory contains the analysis of the Monte_Carlo simulations for the importation time with the two models.
  • analysis_correlation_FOI_diff.py : Create the heatmaps of the importation time for italian provinces.
  • analysis_MC_odd_ratio.py : Compute the odd ratios between diffusion and force-of-infection model.

variables
var  # "d0=0.150R0=0.98" OR "d0=0.150R0=0.98_airport", if inclusion of infection reservoir
which_model  # allfield (force-of-infection) OR diffusion, selection of model
type_net  # regions OR provinces, selection of mobiltiy network
name_parameter # For monocentric: "parameterMC_simu_CHIKV_prov_{name_city}.py", for polycentric: "parameterMC_simu_CHIKV_prov_MULTI.py"

---

Citation

If you use this repository in your research, please cite the associated publication:

@article{YOUR_CITATION_HERE,
  title={},
  author={},
  journal={},
  year={}
}

License

MIT License © 2026 Charley Presigny

Contact

The code was developed and refactored for packaging by Charley Presigny. For questions, bug reports, or collaboration inquiries, contact the corresponding author of the associated publication.