We first summarize the core objects and quantities used throughout the paper.

Our simulations combine a constrained, scale-free–like contact backbone with a multilayer structure and SEIR dynamics. We use a synthetic population (refer to Data and Synthetic Population subsection) to construct a multilayer potential contact network G* and define the H/M groups (refer to Multilayer Potential Contact Network and H/M Classification subsection). We then generate daily contact graphs Gt from G* under the W–m0 framework (refer to Daily Contact Sampling Under the W–m0 Framework subsection) and simulate a discrete-time SEIR model with vaccination and mild or severe branches on these daily networks (refer to SEIR Dynamics With Vaccination subsection). Survey responses are used to calibrate baseline per-layer contact caps (refer to Survey-Based Calibration of Baseline Caps subsection). We report epidemic trajectories and dominance-based outcomes, including the dominance-flip threshold W* (refer to Outcome Measures and Analysis Plan, Baseline Epidemic Dynamics, Dominance Flip Between H and M groups, and Robustness Across Replicate Runs subsections).

Following Ohsawa and Tsubokura [17], we adapt a scale-free network with selfish spatiotemporal constraints as the backbone. In that construction, each node is subject to 2 constraints—an upper bound W on how many others it can meet and a lower cap m0 on how many others it chooses to contact. In our multilayer implementation, these constraints are applied within each layer to capture distinct contact-generating mechanisms. Starting from multiple dense groups, nodes attach new edges while respecting these constraints, producing networks with pronounced communities and a limited number of nodes that bridge communities. In our setting, the nodes are individual agents in a nationally constructed synthetic population, and we embed this backbone into 10 layers (household, school, workplace, distance-driven, and random) that reflect the recorded attributes (refer to Data and Synthetic Population subsection).

Within this backbone, we adopt the W–m0 framework from the study by Ohsawa and Tsubokura [17] to summarize daily contact opportunities at the layer level. W is interpreted as a per-layer daily opportunity upper bound on distinct contacts, and m0 per-layer as a routine-contact minimum on deliberately maintained contacts. We further introduce a random-contact parameter Wrandom that controls the expected number of additional random encounters per person by setting the Erdős–Rényi edge probability on the random layer so that the expected random-layer degree in the potential network G* is approximately Wrandom. In the experiments, we keep per-layer W and m0 fixed and vary Wrandom to represent different levels of random-contact intensity.

A central question in this paper is not only how many infections occur, but which group of individuals contributes most to the spread. To this end, for the high- and medium-contact groups H and M defined on G∗, we track two components: (1) the share of infectious individuals belonging to each group and (2) the average number of effective transmissions generated per infectious individual in that group. Their product denoted ei(W) in the Notation and Key Quantities subsection is served as a dominance score for group i∈{H,M}. As the random-contact parameter Wrandom increases, the dominance scores for H and M can cross, yielding a dominance flip at a threshold W* where eH(W*)=eM(W*). To quantify stochastic uncertainty under fixed parameter settings, we also conduct replicate simulations by varying the random seed and summarize the resulting variability of the estimated dominance-flip threshold W* in the Robustness Across Replicate Runs subsection. The subsequent sections describe how we instantiate the data, networks, sampling scheme, and SEIR dynamics to study this flip.

We use a nationally constructed synthetic population and randomly sample 5000 households (10,038 individuals) from the Tokyo metropolitan area. The synthetic population is generated using a simulated-annealing–based synthetic reconstruction method without individual samples, which algorithmically assigns demographic and household attributes (eg, age, sex, and within-household kinship) to match available population statistics rather than representing real, identifiable individuals [19]. Although this subsample is small relative to the real population, such synthetic populations provide a standard way to reconstruct individual-level contact opportunities while preserving confidentiality [9,20].

Each individual (agent) carries attributes that enable the reconstruction of urban contact layers. An example schema of these individual-level fields is shown in Figure 2. These attributes include:

These attributes support the construction of a multilayer potential contact network on the common vertex set, as described in the Multilayer Potential Contact Network and H/M Classification subsection. In particular, household, school, and workplace layers, together with age-stratified mixing, are chosen to be consistent with empirical contact matrices and their cross-country projections [21,22].

The potential network G* encodes where contacts can occur under the per-layer caps (W,m0,Wrandom) defined in the Notation and Key Quantities subsection. On each day t, we generate a daily contact graph Gt=(V,Et) by sampling a subset of edges Et⊆E from G*. Sampling is carried out for each layer using layer-specific probabilities (Table 2), and all layers are sampled within the same daily step.

For each layer l∈L, each potential edge in that layer is active on day t with probability pl and inactive with probability 1-pl, independently across edges and days. The 4 groups of layers in Table 1 share the same sampling mechanism but use different probabilities pl.

Strongly coupled layers (A)—household, school_classmate, industry_colleague, and best_friends—represent stable, repeated contacts. Their sampling probabilities are set high (Table 2), so that most potential contacts in these layers appear on a typical day.

Organizational layers (B)—industry and school—represent routine organizational contacts such as work and class attendance. Edges in these layers are constructed as maintainable contacts to meet the per-layer routine-contact minimum m0 in the potential network G* and are then activated independently day by day with probabilities pl=0.20. Distance-driven layers (C)—dining_out, amusement, and routine—are constructed with edge probabilities that decay with geographic distance [16,23] and sampled with pl=0.10. Finally, random-layer edges (D) represent chance encounters. Their potential structure is determined by the random-contact parameter Wrandom when G* is built, and on each day, we activate each potential random edge independently with probability pl=0.10. Therefore, the expected realized random-layer degree per day is approximately P_random×W_random. Because this sampling is repeated independently for each day, the sequence of daily graphs {Gt} exhibits natural day-to-day variability even when the underlying potential network G* is fixed.

Using the methodology in the Methods section and the simulation setup and calibration described in the Simulation Setup and Calibration subsection, we examine how changes in contact caps affect epidemic dynamics and the dominance measure ei(W) for the high- and medium-contact groups, H and M defined in the Notation and Key Quantities subsection. The Baseline Epidemic Dynamics subsection presents baseline epidemic curves under the calibrated contact constraints. The Dominance Flip Between H and M Groups subsection then analyzes how the dominance scores eH(W) and eM(W) change as random-contact intensity varies and identifies a dominance-flip threshold. The Robustness Across Replicate Runs subsection evaluates robustness across replicate runs, and the Structural Interpretation of the H→M Dominance Flip subsection interprets this flip in terms of network structure on the multilayer backbone.

Under the baseline random-contact setting (Wrandom=5), the simulated epidemic curves exhibited typical SEIR-type dynamics (Figure 4). Higher overall contact opportunities (W×2, Attr0) produced an earlier and higher infection peak as well as a larger cumulative number of deaths, whereas the half-intensity condition (W/2, Attr2) resulted in slower and more prolonged spread. A more detailed interpretation of how overall contact opportunities affect epidemic speed and size is deferred to the Discussion section.

We next examine how random-contact intensity affects which group of individuals dominates the spread. Throughout this subsection, all layer parameters are fixed to the baseline (Attr1) values in Table 4, and we vary only the random-contact intensity Wrandom while holding the per-layer opportunity upper bound W and the per-layer routine-contact minimum m0 fixed (for the nonrandom layers).

As defined in the Notation and Key Quantities subsection, the dominance score for group i∈{H,M} at a given contact-opportunity setting W is

where si(W) is the share of infectious individuals belonging to group i and ni(W) is the average number of effective transmissions generated by one infectious individual in that group. In the simulations below, we evaluate ei(W) at different values of the random-contact intensity Wrandom, using W as the contact-opportunity argument.

For each tested value of Wrandom between 5 and 10, we run multiple epidemics under the baseline (Attr1) setting and compute the dominance scores eH(W) and eM(W) from the simulated trajectories. We then fit smooth curves to their dependence on random-contact intensity. A dominance flip occurs when the dominant group changes between H and M as W varies (refer to Notation and Key Quantities subsection). Under the baseline condition, the fitted curves cross near W*≈9.82.

As we will next show in the Robustness Across Replicate Runs subsection, repeated runs yield a narrow but nonzero range for this estimate (median 9.8294, IQR 9.8231‐9.8367), so we report W* as an estimated range within this setting rather than a fixed constant. For example, when Wrandom=9, the dominance scores are approximately eH≈0.37 and eM≈0.30, so the high-contact group H still dominates. When Wrandom=10, they become eH≈0.32 and eM≈0.33, so the medium-contact group M becomes dominant. These values illustrate the transition from an H-dominant regime to an M-dominant regime as random-contact intensity increases.

Figure 5 summarizes this pattern through the difference eH(W)-eM(W). The horizontal axis is parameterized by the random-contact intensity Wrandom, and the sign of eH(W)-eM(W) changes at the threshold W*≈9.82; it is positive for Wrandom≲9.8 (H-dominant) and negative for Wrandom≳9.8 (M-dominant).

The daily contact sampling and epidemic progression in our simulations are stochastic. To assess how this stochasticity affects the estimated dominance-flip threshold, we conducted replicate simulations under identical model parameters while varying the random seed. In each replicate, we estimated the dominance-flip threshold W* using the same definition and estimation procedure as in the Notation and Key Quantities and Dominance Flip Between H and M Groups subsections, based on the crossing point of the dominance scores ei(W)=si(W)ni(W).

Across replicates, the estimated W* values were tightly concentrated (Figure 6). The median estimated threshold was 9.8294 (IQR 9.8231-9.8367). Using the same tested grid of the random-contact intensity parameter (Wrandom∈ [5,10]), a dominance-score crossing within this range was observed in 30 out of 30 replicates. These results indicate that the dominance flip is consistently observed under the fixed parameter setting, and that the estimated threshold has a measurable uncertainty range.

The previous subsection showed that, under the W–m0 framework, the dominance score ei(W) for the high-contact group H and the medium-contact group M crosses at a threshold W*≈9.8 in terms of the random-contact parameter Wrandom. We now examine how the multilayer structure of the network changes around this threshold. We use a same-node visualization of a representative subgraph to compare the allocation of intercommunity bridges before and after the flip. This same-node visualization is intended as a mechanistic illustration rather than as an estimator of W*; robustness of the flip threshold is assessed separately through replicate runs.

From the realized daily contact graphs on a fixed simulation day, we extract a same-node subgraph (Figure 1) anchored on high-contact individuals. We first identify the 3 highest-degree nodes in group H, take their 1‐2 hop neighborhoods, and cap the size at about 250 nodes. To provide a stable reference for both panels, we detect communities (modules) on the aggregated edge union graph of the 2 panels (same node set) using a standard label-propagation algorithm, yielding a nonoverlapping partition that is used only for visualization and for defining cross-community edges. An intercommunity bridge is an edge whose end points belong to 2 different communities. Among these cross-community edges, we highlight a bridge backbone ranked by edge betweenness centrality; higher-ranked edges are drawn thicker and in orange. Red nodes indicate high-contact individuals (H), and blue nodes indicate medium-contact individuals (M).

The left panel in Figure 1 corresponds to a setting with weaker random contact (Wrandom§amp;lt;W*), where H dominates. Here, most of the thick orange bridges connect adjacent background regions and have at least 1 red end point; many visually prominent bridges start from the central red core and end at red nodes in nearby communities. Blue nodes appear on bridges mainly as peripheral end points at the edges of communities. The right panel corresponds to stronger random contact (Wrandom§amp;gt;W*), where M dominates. In this case, orange bridges are more numerous and many of the thick orange edges run across nonadjacent background regions, linking visibly distant parts of the layout. Compared with the left panel, a larger fraction of these thick bridges now have blue end points, and several orange edges connect 2 blue nodes. Thus, medium-contact individuals more often sit at the ends of intercommunity bridges and connect a larger number of distinct communities. This reallocation of bridges increases the effective reach of the groups M and raises eM(W) relative to eH(W), consistent with the dominance flip observed in the Dominance Flip Between H and M Groups subsection.

Within the multilayer, daily-resampled contact network and discrete-time SEIR model described in the Methods and Results sections, we examined how overall contact opportunities and random-contact intensity shape epidemic outcomes and the leading contributors to spread. First, under fixed random-contact intensity (Wrandom=5), increasing the opportunity upper bound W (twofold vs baseline vs half; Figure 4) produced earlier and higher peaks and larger cumulative deaths. Thus, within the calibrated range of parameters, higher overall contact opportunities directly increased both the speed and scale of epidemic expansion. Second, focusing on the standard overall-intensity condition and varying the random-contact parameter Wrandom, we observed a dominance flip between high-contact (H) and medium-contact (M) individuals: for low Wrandom the H group dominates, whereas for high Wrandom the M group dominates (refer to Dominance Flip Between H and M Groups subsection; Figure 5). Dominance is measured by the score ei(W)=si(W)ni(W) defined in the Notation and Key Quantities subsection, and the empirical threshold at which eH(W)=eM(W) occurs is Wrandom*≈9.8. These results show that changes in random-contact intensity can reshuffle which part of the population contributes most to transmission, even when overall contact opportunities are held fixed.

The dominance flip arises from the interaction of 2 effects, which we summarize here and illustrate using 2 complementary views (Figures 1 and 7).

To provide an intuitive interpretation within the same contact-constraint framework, Figure 7 illustrates a schematic to meet the routine-contact minimums (W,m0). Individuals are ordered on the horizontal axis by baseline contact intensity on the potential network G* (left: higher-degree nodes, mainly H; right: lower-degree nodes, mainly M), and the vertical axis represents expected daily contacts. As random-contact intensity Wrandom increases while W and m0 are held fixed, the curves shift upward. Following Ohsawa and Tsubokura [17], the dashed lines mark m0 and 2m0 as reference levels. In this schematic, increasing Wrandom raises contacts more on the M side where residual opportunity under (W,m0) remains, so more M nodes cross the 2m0 level, while additional contacts on the H side are limited by saturation. This qualitative pattern is consistent with the observed increase in nM(W) relative to nH(W) and the resulting dominance flip reported in the Dominance Flip Between H and M Groups subsection.

This mechanism is consistent with constrained-interaction views that separate opportunity (W) from routine-contact minimum (m0) in the W–m0 framework [17], and with previous work showing that short-time reshuffling and intercommunity bridges strongly influence invasion routes and timing in modular networks [10,11].

Methodologically, this study combines 3 components in a single, interpretable framework. First, the W–m0 setup separates a daily opportunity upper bound W from a routine-contact minimum m0, and introduces Wrandom as a separate control for random encounters (refer to Notation and Key Quantities subsection). Second, these quantities are implemented on a multilayer contact backbone with community structure and realistic venue types (household, school, workplace, distance-driven layers, and a random layer; refer to Model Overview, Data and Synthetic Population, Multilayer Potential Contact Network and H/M Classification, and Daily Contact Sampling Under the W–m0 Framework subsections). Third, an extended SEIR model with vaccination and mild or severe branches is run on the daily graphs {Gt}, with contact-related parameters calibrated from an online survey in Tokyo and Kanagawa (refer to Survey-Based Calibration of Baseline Caps subsection). Within this setting, a small set of dominance-related quantities, {si(W),ni(W),ei(W),W*}, summarizes which degree-based group (H or M) contributes most to transmission under each contact scenario. Thus, the framework retains structural richness while keeping the key indicators of “who leads spread” simple and comparable across settings.

Substantively, the results refine common intuitions about targeting high-contact individuals. Many models and narratives implicitly assume that high-contact individuals are always the main drivers of spread. In our multilayer, community-structured network, this holds when random-contact intensity is low, for Wrandom§amp;lt;Wrandom* the dominance score is larger for H. However, as random-contact intensity increases and Wrandom passes the empirical threshold Wrandom*≈9.8, the medium-contact group M becomes dominant. This shows that changes in random-contact opportunities can flip which part of the population contributes most to transmission, even when the overall opportunity bound W and the routine-contact minimum m0 are fixed. This result complements previous bridge-focused and high-degree targeting narratives by showing that who occupies effective bridging positions can change when incidental encounters increase, even if per-contact transmissibility is unchanged.

These findings have direct implications for intervention priorities. When Wrandom§amp;lt;Wrandom*, most effective transmissions originate from the H group. In this regime, measures that concentrate on H-type settings—such as large workplaces, dense households, and other high-degree environments—are expected to yield the largest marginal reductions in transmission (eg, focused screening or vaccination and testing outreach around identified H clusters). As Wrandom approaches or exceeds Wrandom*, random encounters play a larger role and the contribution of M increases. In this regime, reducing random-contact opportunities becomes important. In our framework, such measures aim to decrease the effective value of Wrandom or of the gap W-m0, for example, by limiting the size and frequency of large events, spreading flows over time through staggered commuting or opening hours, and reducing dense queues or long stays in crowded transfer spaces. Viewed through the dominance metric, these measures limit the formation of new intercommunity bridges and help prevent a shift from H-dominant to M-dominant spread.

Several limitations should be noted.

First, the synthetic population is restricted to a subsample of the Tokyo metropolitan area, and the layer-wise parameters were calibrated using 1 online survey in Tokyo and Kanagawa (refer to Survey-Based Calibration of Baseline Caps subsection), which is generally limited in terms of generalizability. This design may not represent typical daily behavior and may limit generalizability. The study is therefore not intended as a full reproduction of real-world contact patterns across settings. Instead, it is a scenario analysis focused on higher-activity conditions, where contacts are elevated and dominance shifts can be examined more clearly. Low-activity periods are outside the scope of this work.

Second, the quantitative value of the threshold Wrandom* depends on epidemiological parameters, such as infection probabilities and progression rates. Our structural conclusion—that increasing Wrandom in a heterogeneous, multilayer network can flip dominance from H to M—relies on the presence of degree and layer heterogeneity together with an increasing Wrandom, rather than on a specific point estimate of any single parameter. We address uncertainty from randomness in daily contact sampling and transmission by reporting replicate runs under fixed parameter settings. However, these replicate runs do not replace a full sensitivity analysis over epidemiological assumptions (eg, infection rates, vaccine effects, and severity parameters).

Third, the operational definitions of the random-contact layer and of the H/M groups may differ across settings. In this study, H and M are defined by degree quantiles on the potential network, and the random layer aggregates several types of chance encounters; other implementations that respect the W–m0 structure are possible.

Finally, our results are averages over a finite number of stochastic realizations. Even under the same model setting, the estimated value of Wrandom* varies across simulation runs due to randomness in daily contact sampling and transmission events. Importantly, while the point estimate of Wrandom* shifts across realizations, the qualitative tendency for an H→M dominance flip as Wrandom increases remains robust.

In this study, we modeled urban daily life as a multilayer contact network with a constrained, scale-free–like backbone and SEIR dynamics with vaccination. Individual contacts were governed by per-layer opportunity caps and per-layer routine-contact minima, together with a separate parameter for random encounters, and these quantities were calibrated using survey data from Tokyo and Kanagawa. On this backbone, we followed the epidemic on daily resampled networks and evaluated which degree-based group—H or M individuals—contributed most to transmission.

Within the range of calibrated settings, higher overall contact opportunities led to earlier and higher epidemic peaks and larger cumulative deaths, while lower opportunities slowed and prolonged spread. Beyond this expected pattern, our main finding is that changing random-contact intensity, even when overall contact opportunities and routine-contact minimum are held fixed, can flip which group dominates the spread. When random encounters are limited, most effective transmissions arise from the high-contact group. As random contact increases, saturation on the H side and the proliferation of intercommunity bridges through medium-contact individuals together shift the main contribution to the M group. In other words, who leads spread depends not only on baseline contact opportunities, but also on how many of those contacts are chance encounters that connect different communities. Replicate simulations further quantify stochastic uncertainty under fixed settings, showing that W* varies within a measurable range while the flip is consistently observed.

These results suggest that intervention priorities should reflect the level of random contact. When random encounters are relatively rare, focusing resources on high-contact settings—such as large workplaces and dense households—is expected to be most effective, because most transmission chains are rooted in these environments. When random encounters are frequent, reducing chance encounters and cross-community mixing becomes more important. In this framework, such measures can be interpreted as lowering the effective intensity of random contact or narrowing the within-layer gap between opportunity caps and routine-contact minima, for example, by limiting crowding and dwell time in shared spaces, staggering flows over time, and avoiding prolonged close contact with people outside one’s usual circles.

This work has several limitations and is one step in a longer line of research. The synthetic population and layer parameters are based on one metropolitan area and a COVID-19–like infection, and the classification into high- and medium-contact groups is only one possible way to summarize heterogeneity. Network dynamics are driven by daily resampling under fixed caps, without explicit behavioral feedback or seasonal change. Future studies should apply the same framework to other regions and time periods, investigate how the dominance flip and its threshold vary under different epidemiological and behavioral conditions, and combine additional data sources to refine the calibration of contact opportunities, routine-contact minimum, and random encounters. Even with these simplifications, the results of this study indicate that distinguishing between routine and random contacts, and tracking which group dominates spread, can provide a simple and interpretable guide for thinking about targeted control in multilayer urban networks.