Multiple contact settings

CountryData.contact_matrix accepts either a single Matrix{Float64} or a Vector{Matrix{Float64}}. Passing a vector lets the model treat distinct social contact regimes. For example, community mixing and workplace mixing — as separate matrices that are combined inside the ODE force-of-infection calculation.

Creating a multi-setting CountryData struct

Retrieve a country struct, then replace its contact_matrix field with a vector of matrices:

using Daedalus

cd = Daedalus.DataLoader.get_country("Australia")
cm = deepcopy(cd.contact_matrix) # baseline 4×4 contact matrix
cm_work = cm .* 2.0 # a second setting with doubled contacts

cd_multi = deepcopy(cd)
cd_multi.contact_matrix = [cm, cm_work]

2-element Vector{Matrix{Float64}}:
 [3.71875 2.5982167689627107 5.739112393931974 0.2728100904551632; 0.909536863622831 13.062337428330006 5.741991759560344 0.5229290973202789; 0.6420045383488484 1.8348941187717769 11.256654747715649 1.0003494672654953; 0.13475822848175384 0.6540519090439044 3.760931041766723 2.542189501604722]
 [7.4375 5.196433537925421 11.478224787863947 0.5456201809103264; 1.819073727245662 26.12467485666001 11.483983519120688 1.0458581946405578; 1.2840090766976968 3.6697882375435538 22.513309495431297 2.0006989345309907; 0.2695164569635077 1.3081038180878088 7.521862083533446 5.084379003209444]

Data.get_settings returns the number of active contact settings:

Daedalus.Data.get_settings(cd)

1

Daedalus.Data.get_settings(cd_multi)

2

Running the model with multiple settings

Pass the modified struct directly to daedalus along with infection parameters:

result_multi = daedalus(cd_multi, "sars-cov-2 delta", time_end = 300.0)

times  = Daedalus.Outputs.get_times(result_multi)
deaths = Daedalus.Outputs.get_values(result_multi, "D", 1)
println("Total deaths at end: ", round(Int, last(deaths)))

Total deaths at end: 190293

Calculating β with multiple settings

daedalus computes β from the sum of all contact matrices via Data.total_contacts:

# NOTE: example, this code is not run
contacts_unscaled = total_contacts(prepare_contacts(cd; scaled = false))
beta = get_beta(contacts_unscaled, r0, ...)

The effective reproduction number is therefore calibrated against the total contact rate across all settings. As a consequence, splitting contacts across multiple equal settings does not change the epidemic. Adding a second identical matrix doubles the total contacts, which halves β.

cd_double = deepcopy(cd)
cd_double.contact_matrix = [cm, cm]

infection = Daedalus.DataLoader.get_pathogen("sars-cov-2 delta");

# get beta used for each country
beta_single = Daedalus.Helpers.get_beta(
    Daedalus.Data.total_contacts(
    Daedalus.Data.prepare_contacts(cd; scaled=false)
), infection)

beta_double = Daedalus.Helpers.get_beta(
    Daedalus.Data.total_contacts(
    Daedalus.Data.prepare_contacts(cd_double; scaled=false)
), infection)

0.0016750668855633447

Print $\beta$ values to check that beta_double is half of beta_single.

println("β (1 setting): $beta_single")
println("β (2 settings): $beta_double")

β (1 setting): 0.0033501337711266893
β (2 settings): 0.0016750668855633447

How settings appear in the ODE

The matrices are stacked into a 3D array of shape (N_TOTAL_GROUPS, N_TOTAL_GROUPS, K) where K is the number of settings. In each ODE evaluation, a weight vector of length K contracts the third dimension via Helpers.weighted_slice_sum!:

c3d = Daedalus.Data.contacts3d(cd_multi)
size(c3d) # (49, 49, 2) for two settings

(49, 49, 2)

Current weights are ones(K) (all settings fully active). This design is intended to support future NPI scaling, where individual settings — for example, workplace contacts during an economic closure — can be reduced independently without affecting community contacts. ```