Accurately forecasting population dynamics, driven by natural growth and mobility, is essential for effective urban planning and territorial development. Traditional methods, such as Cohort Component Analysis (CCA), primarily rely on historical data and often overlook the nonlinear interdependencies among multiple demographic parameters and interacting regions. Consequently, these rigid frameworks can lead to outdated predictions and overlooked implications (Rees and Wilson, 1977). To overcome these limitations, this study benefits from the dynamic modeling capabilities of an enhanced Lotka-Volterra model for human population prediction, while simultaneously integrating a refined probabilistic gravity approach. This integrated framework offers a more precise and comprehensive tool for understanding and predicting population dynamics, thereby improving the accuracy of regional planning and policymaking (Lutz and Gailey, 2020).

Population estimates generally rely on the evolution of various demographic components and prevailing trends (Reinhold and Thomsen, 2015). CCA, pioneered by Whelpton (1928) in Canada, remains widely employed for population enumeration and projections. It is widely adopted in developed nations for its effectiveness in capturing historical trends and the evolution of natural parameters such as birth and death rates, classified by age and gender (Chi, 2009; Wilson, 2016). CCA employs straightforward methodologies, assuming linearly independent impacts of typical demographic attributes while often underrepresenting population mobility and other abstract parameters influencing regional dynamics (Smith et al., 2013). The complexity of demographic forecasting essentially lies in interpreting these interconnected mobility patterns (Lucas, 2015; Lutz and Gailey, 2019; Stawarz et al., 2021; Murayama and Nagayasu, 2021).

Interregional mobility is significantly influenced by individuals’ pursuit of improved housing, job opportunities, vibrant communities, and higher social and environmental standards (Greenwood, 1997; Teixeira, 2009; Seo et al., 2020). Anticipating its socioeconomic consequences and demographic shifts is essential for effective planning, management, and territorial development (Champion and Fielding, 1992; Zoraghein and O’Neil, 2020). Critical decisions on social welfare, economic resource allocation, land use planning, property valuation, taxation, and public services, especially transportation, depend on robust demographic insights (Biswas et al., 2009; Terroso-Sáenz and Muñoz, 2022). However, a static and reactive approach in isolated regions fails to harness demographic trends effectively, exacerbating issues that require strategic intervention for balanced regional development (Lichter and Johnson, 2006). Such approaches often overlook dynamic socioeconomic factors, leading to fragmented and suboptimal development outcomes (Miao and Pojani, 2022). Integrated and adaptive planning methods are essential for responding to changing environmental and economic conditions, as rigid, isolated strategies can hinder effective regional development (Özdilek, 2016; Wang and Li, 2021).

Demography and migration inherently encompass socioeconomic and geographic dimensions, making them amenable to modeling cultural dynamics (Voss, 2007; Burton et al., 2010; Matthews and Parker, 2013; Wilson et al., 2021). To refine population growth estimations, demographers have traditionally utilized age-sex spatio-temporal life tables or actuarial probability matrices, sometimes integrating origin-destination mobility data (Rau et al., 2022; Malmberg et al., 2023). While many studies focus on internal migration rates, interactions between neighboring regions—driven by geographical disparities and socioeconomic ties—catalyze complex inter-regional dynamics (Song et al., 2013; Huang et al., 2018). The integration of spatial interaction models, particularly those utilizing entropy-maximizing principles as established by Wilson (1970), Wilson (2010), provides a robust framework for incorporating spatial dependencies of socioeconomic and migration patterns.

In recent years, advanced computational algorithms have significantly enhanced demographic mobility estimates by integrating interregional dependencies. Techniques such as feedforward neural networks (Openshaw, 1993), agent-based modeling (Billari and Prskawetz, 2003; Gopal, 2006), the radiation model (Simini et al., 2021), and Graph Convolutional Network (GCN) architectures (Wang et al., 2023a) have expanded the analytical toolkit. These methodologies, ranging from time-series analysis to Bayesian modeling and machine learning approaches, improve the predictive capabilities of demographic studies (Abel et al., 2013; Robinson and Dilkina, 2018; Sohst and Tjaden, 2020). By accounting for spatial and temporal interdependencies and the complex interplay of socioeconomic factors, these advancements offer a more nuanced understanding of population dynamics.

Various tools and specific programs enhance predictive capabilities by integrating spatial and temporal data, making them indispensable for modern demographic research and planning (Wang et al., 2023b). These tools leverage programming and simulation approaches to anticipate territorial development, creating virtual spaces that facilitate hypothesis testing and scenario development (Deng et al., 2023). For instance, Origin-Destination REIMHS analysis supports traffic flow and population movement predictions, providing insights into regional interactions (Wang et al., 2019). MEPLAN, a land-use and transport interaction model, integrates the economic impacts of transportation systems, offering a holistic view of regional development (Pignataro et al., 1998; Hunt et al., 2005; Cooke and Leydesdorff, 2006). The PECAS model further refines understanding by simulating economic and spatial interactions in urban areas, improving population predictions (Ellam et al., 2018). Additionally, TELUM and UrbanSim models simulate household and workforce location behaviors under various land use scenarios, aiding in the analysis of demographic trends and their impact on regional planning (Waddell, 2002; Grant, 2002; Deng et al., 2023).

The global evolution of biological systems and other complex phenomena necessitates nonlinear dynamical models to capture space-time variations (Hofbauer and Sigmund, 1998; Burch, 2018). The Lotka-Volterra model, employing differential equations within a prey-predator dynamic framework, has found applications in epidemiology and economics (Cockshott and Cottrell, 2009; Hol et al., 2016). However, its direct use in human population prediction remains less explored (Bischi et al., 2000). Despite its theoretical robustness, demographic forecasting often relies on accounting methods lacking a solid theoretical basis (Hudson, 1970; Tranos et al., 2015a; Tranos et al., 2015b). The Lotka-Volterra model’s mathematical foundation offers a promising avenue for advancing demographic forecasting, providing a coherent framework for understanding population dynamics and guiding effective resource allocation (de Valk et al., 2022). Although not originally a demographic model, the gravity model, based on Newton’s law of gravitation, has been applied to study migration between competing regions (Reilly, 1931; Zipf, 1949; Anderson, 2011), but requires adjustments for social science contexts (Davis et al., 2013).

Integrating the Lotka-Volterra dynamic equations enables models to more effectively capture interactions between competing populations. Additionally, advanced gravity models, particularly those based on entropy-maximizing principles (Wilson, 2008), can enhance population prediction capabilities. However, even integrated Lotka-Volterra and gravity models do not fully account for all the complexities of population mobility and predictability. Simplifications of parameters and phenomena make accurate predictions of social behaviors, migration patterns, and population growth challenging (Partridge et al., 2007). Moreover, factors such as subjective decision-making, cultural influences, and policy changes further limit the comprehensive predictive power of these approaches (Billari and Prskawetz, 2003).

This study presents an innovative methodological framework to address the intricacies of population dynamics across three regional clusters in Quebec, Canada. While it focuses on regions within the Quebec province, the presented model—its concepts, data, and methodology—is sufficiently generalizable to other geographic contexts, provided suitable data and identified clusters. The integrated model’s ability to capture dynamic interactions makes it highly adaptable. Future studies could apply this approach elsewhere, considering local data and specific regional dynamics.

By incorporating concepts such as carrying capacity, multiple interacting regions, and resource constraints, and by fine-tuning parameters using historical data, our adapted model closely simulates competitive and cooperative dynamics reflective of real-world interactions. This comprehensive approach enhances the accuracy and relevance of interregional migration fluidity and provides a robust tool for demographic forecasting. Demonstrated through illustrative examples and optimized using advanced computational techniques, our model offers reliable, policy-informing insights into regional population influences over time, significantly advancing demographic forecasting and supporting strategic planning and effective policy formulation.

Effective demographic predictions hinge on concepts, data, and mathematical approaches grounded in fundamental theories (Bryant and Zhang, 2019). Current data-centric methodologies in demography have not significantly advanced the field (Luhmann, 1995; Burch, 2018), largely due to the absence of a comprehensive foundational framework. This often leads to the ecological fallacy, where aggregate data is misapplied to individual behavior, failing to reflect real-world complexities (Courgeau, 2007; Courgeau, 2012). Consequently, the field has been slow to develop robust theories that integrate social, economic, and geographic variables, frequently drawing on insights from other disciplines without forging a cohesive theoretical framework (Doran, 2000; Silverman et al., 2011; Silverman et al., 2013; Courgeau et al., 2017a). This lack of a unified approach hampers the ability to fully capture the nuances of population dynamics (Burch, 2002).

The literature underscores the significance of concepts and demographic parameters as the foundation for various methods used to study population (Paul and Yeoh, 2020). Parameters such as aging and birth/death rates, governed by biological rules, are relatively predictable (Guillot, 2011). However, parameters related to human judgment and behavior—such as fertility, lifestyle, family composition, migration, and regional socioeconomic disparities—are less predictable yet crucial (Kuate Defo, 2014). For instance, individuals often relocate from high-unemployment areas to regions offering lower unemployment rates (Myrdal, 1957; Polese, 1981). Lee’s migration theory (1966, 1969) elucidates push and pull factors by considering both attractive and repulsive regional attributes. These components span biology, sociology, economics, geography, politics, and physics (Rogers and Castro, 1981; Rogers and Willekens, 1986; de Valk et al., 2022).

The economic impacts of utility parameters, explained by Lancaster (1966) and quantified by Rosen (1974), provide a framework for understanding decision-making through concepts, data, and methods. Just as rational economic agents drive markets toward equilibrium, individuals make demographic decisions that influence population variations via migration dynamics (Bloom et al., 2009; Aksoy et al., 2019). For instance, unemployment significantly affects migration patterns (Wiginton, 2013). Tiebout’s (1956) theory suggests that “people vote with their feet,” indicating preferences for locations with better job availability and conditions (McDonald, 2004; Walsh et al., 2011). Cultural and social disparities highlight mobility dynamics, particularly from poorer to wealthier areas (Kritz et al., 2013), as originally noted by the Chicago School (Harris and Ullman, 1945). City size plays a crucial role, with larger cities typically attracting more immigrants due to broader opportunities and amenities, while smaller cities may attract fewer migrants due to limited resources (Malecki, 2004; Hyndman et al., 2006; Price and Benton-Short, 2008). Once settled, immigrants often move from their initial cities to larger urban centers (Newbold, 2007; Newbold, 2011). Governments engage in optimally mitigating immigration policies to balance total population evolution and ensure service provision and economic growth (Green and Green, 2004).

Geographers and urban planners apply geographical theories to assess location utility and rent (Alonso, 1960; Oates, 1969). Tobler’s (1970) First Law of Geography states that “everything is related to everything else, but near things are more related than distant things.” Spatial proximity and interregional relations affect flow intensity (Ullman, 1954), typically measured by geographical distance (Fotheringham, 1981). Migration flows represent geographically embedded networks whose intensity correlates with population movement (Lee, 1966; 1969; Batten and Boyce, 1987; Abel and Sander, 2014). Demographers use a spatio-temporal probability matrix to integrate mobility resulting from observed trends. The stationary regional populations of such a life table serve as inputs to numerical demographic calculations. Based on various age-structured data like birth location and actual residence, these matrices improve population growth estimates (Leslie, 1945; Rogers and Castro, 1981). In this perspective, geographers and demographers draw on Newton’s law of gravitation to explain interaction and migration over space (Reilly, 1931; Zipf, 1949). The most commonly applied form of the gravity model predicts migration interaction between two locations of differing importance, factoring in distances between them (Poot et al., 2016).

The gravity model aligns well with migration theories and offers straightforward predictive capabilities for population growth (Stewart, 1950; Anderson, 2011). Based on Newton (1687) Law of Universal Gravitation, the model suggests that regional attractiveness increases flow intensity, with distance acting as a cost factor (Haynes and Fotheringham, 1984). Enhancements to the gravity model incorporate additional regional attributes (Lewer and Van den Berg, 2008) and advanced iterations that account for multiple push/pull factors and spatial spillovers (Fik and Mulligan, 1998). These refinements consider internal, interstate, and international migrations, providing a comprehensive analytical tool for demographic studies (Karemera et al., 2000; Peeters, 2012; Beine et al., 2016).

Recent advancements in demographic research underscore the importance of integrating dynamic models and computational techniques to better understand and predict population behaviors. Burch (2018) emphasizes moving beyond traditional data-centric approaches, advocating for models that capture the complexity of human decision-making and interactions. Agent-based models and system dynamics serve as essential tools for integrating micro- and macro-demographic phenomena, offering a comprehensive framework for analysis (Billari and Prskawetz, 2003). Courgeau et al. (2017b) suggest that adopting these methods can lead to broader, model-based research programs with enhanced explanatory power. The use of mathematical modeling and simulation tools, such as the Lotka-Volterra model, is crucial for reflecting non-linear interactions and dynamic mobility within populations (Franck, 2002).

Dynamic models describe the temporal evolution of physical, ecological, economic, financial, and social systems (Puu, 2000; Bischi et al., 2000). Discrete time dynamic models, such as the logistic map, exhibit complex, chaotic behavior characterized by trajectories converging to attractor points or basins of attraction (Gumowski and Mira, 1980; Bischi and Tramontana, 2010; Özdilek, 2022). The coexistence of multiple attractors necessitates spatial delineation and equilibrium analysis, illustrating path-dependent processes where long-term outcomes depend on initial conditions (Liu and Xiao, 2006; Liu and Xiao, 2007; Blackmore et al., 2001). Such complexities necessitate comprehensive global dynamical modeling analyses, exemplified by the Lotka-Volterra equations, akin to those used by Lorenz (1963) in weather prediction.

Smale (1976) demonstrated how ordinary differential equations describe the dynamics of competing species, supported by Rescigno and Richardson (1967) and May (1973). Gilpin and Feldman (2017) showed how prey evolution in response to predation affects population dynamics, influenced by temporal fitness landscapes. Although natural selection progresses more slowly than population dynamics, demographic shifts triggered by predator invasions or environmental changes can occur swiftly (Hansen et al., 2013; Cortez and Weitz, 2014). Thus, classical prey-predator models must incorporate eco-evolutionary processes (Mougi and Iwasa, 2010). Even simple deterministic two-species models can exhibit unstable cycling and chaotic interactions, as seen in microbial communities following perturbations like pesticide treatment or pathogen colonization (Hol et al., 2016; Remien et al., 2021).

The Lotka-Volterra framework posits that growth rates depend on the quantities of other components. In a fox-rabbit system, for example, the prey consumption rate by predators correlates with prey abundance, influencing fox populations and vice versa. This model assumes a fixed territory and exclusive prey, theoretically allowing indefinite predator growth as prey decline, highlighting an inherent power imbalance (Burch, 2018). Beyond ecology, prey-predator models have explained competitive, cooperative, and predatory strategies in various contexts, including industry growth patterns and resource competition (Bischi and Tramontana, 2010; Dendrinos and Mullally, 1985; Nijkamp and Reggiani, 1992; 1998; Pijush et al., 2021).

The original Lotka-Volterra equations, describing predator-prey dynamics through nonlinear first-order ordinary differential equations and parameters like birth and death rates, can be adapted to complex real-world relationships. By modeling populations competing for common resources such as education, jobs, and housing, these equations can extend to human demographic projections (Lotka, 1925; Volterra, 1926). Although initially rooted in animal population dynamics, careful justification and adaptation enable examining human interactions and resource competition, offering insights into population trends and dynamics.

While few studies have applied Lotka-Volterra models to human populations and socioeconomic contexts, most have not integrated spatial mobility or focused on detailed population forecasting. For instance, Bischi and Lamantia (2002) modeled competition and cooperation among agents in natural resource exploitation without emphasizing spatial aspects. Izquierdo and Hannachi (2018) used Lotka-Volterra equations for urban dynamics, emphasizing complex systems theory but not integrating gravity models or predictive demographic modeling.

The gravity model, extensively employed in migration studies, represents spatial interactions based on population sizes and distances between regions (Wilson, 1970). Although Batty (2005) and Stillwell et al. (2014) highlighted its role in understanding urban growth and internal migration, these applications generally operate independently of dynamic population models. Attempts to combine Lotka-Volterra and spatial interaction models have provided theoretical insights but not a unified framework for practical human population prediction (Wilson, 2006; Wilson, 2008; LeSage and Pace, 2009; Wilson and Dearden, 2011). As a result, existing models, while insightful, do not fully address complexities like population mobility and predictability. Simplified parameters, subjective decision-making, cultural influences, and policy changes remain challenging to incorporate, limiting comprehensive predictive power (Billari and Prskawetz, 2003).

Recent contributions emphasize that fully probabilistic and Bayesian frameworks can rigorously capture uncertainty, moving beyond deterministic trajectories toward probability distributions of future demographic indicators (Keilman, 2018). Such forecasts, as illustrated by Raftery and Ševčíková (2023) and adopted by international agencies, provide predictive intervals clarifying likely outcomes. Similarly, Wiśniowski et al. (2015) enhance the Lee-Carter model with a Bayesian hierarchical approach, generating posterior distributions rather than single estimates. At smaller spatial scales, Wilson et al. (2021) highlight the need for probabilistic techniques to handle volatility and complexity, while Azose et al. (2016) incorporate migration uncertainty so that all major demographic components—fertility, mortality, and migration—are probabilistically modeled.

These Bayesian and probabilistic frameworks excel at uncertainty quantification, relying on statistical and data-driven foundations. By contrast, this study employs a theory-driven, dynamic model integrating Lotka-Volterra equations with a probabilistic gravity concept to capture nonlinear interdependencies and spatial interactions. While the current approach produces interpretable trajectories and equilibrium points, it has yet to provide full probability distributions. Future work could incorporate Bayesian inference for parameter estimation and credible intervals, thereby melding the strong theoretical underpinnings with robust uncertainty quantification, creating a more comprehensive and adaptive forecasting tool.

The population projection widely uses the Cohort Component Model (CCM). Statistics Canada commonly predicts populations 1 year at a time by extending past trends based on CCM, with slight differences between provinces. The CCM adopts a top-down population estimation at the national level first, followed by subnational projections, producing projections by age and sex (Cameron et al., 2007; Cameron and Poot, 2010; Cameron and Poot, 2011). The CCM’s Equation 1 rests on objective counting of individuals in a population to sequentially build two basic components in recursive form: P t + i = P 0 + P BD + P IO (1) where P t + i is the population estimate at time t + i ; P 0 is the initial population, P DB is the natural increase or the balance, for instance between birth and death of population, and P IO is the net internal in/out-migration or mobility which might include interprovincial and international migration. The first term P BD ​ is quite well predictable with several demographic parameters like the future fertility and mortality rates, benefiting from various modeling approaches improving the accuracy of each of these parameters (Archambault, 1999; Lapointe et al., 2006). However, the precision of population projections mainly relies on the second component P IO , which is rather complex and less predictable (Dunn, 2006).

Traditional demographic methods, such as the CCM, often fail to capture the nonlinear interactions and dynamic spatial dependencies inherent in population movements (Rees and Wilson, 1977; Wilson, 2010). To overcome these limitations, this study employs the Lotka-Volterra prey-predator model, which effectively models competitive and cooperative interactions among populations. The original equations of the Lotka-Volterra model are: dx t dt = ax t − g 1 x t y t (2) dy t dt = − cy t + g 2 x t y t (3) where x t represents the resource (prey) population and y t represents the consumer (predator) population. In these equations, a is the natural exponential growth rate of the prey x in the absence of predators, indicating that the prey population grows exponentially when not constrained by predation. The term g 1 x t y t in Equation 2 represents the rate at which the prey population is consumed by predators, with g 1 being the interaction coefficient that quantifies this effect.

The differential Equation 3 describes the dynamics of the predator population. Here, c represents the natural mortality rate of the predators in the absence of prey, indicating a decline in the predator population over time without a food source. The term g 2 x t y t captures the growth of the predator population due to predation, with g 2 being the interaction coefficient that quantifies how effectively the predator population increases based on the availability of prey.

In the absence of a predator ( g 1 = 0 ), Equation 2 simplifies to: dx t dt = ax t and dx t x t = adt

Integrating this differential equation, we get the solution: x = K e at (4)

This Equation 4 describes the exponential growth rate of prey in the absence of predators, leading to an indefinite increase in the prey population (birth process).

For the rate of decline of predators, Equation 3 in the absence of prey ( g 2 = 0 ) simplifies to dy t dt = − cy t and dy t y t = − cdt . Integrating this differential equation, we obtain: y = Ke − ct (5)

This Equation 5 describes the exponential decline of predators in the absence of prey (mortality process).

To adapt the Lotka-Volterra model for human demographic applications, we incorporate a carrying capacity condition, as originally proposed by Verhulst (1838), to account for environmental constraints on population growth. As shown in Equation 6, the logistic growth model describes how populations expand in environments with limited resources by incorporating the carrying capacity K, which represents the maximum population size an environment can sustain indefinitely. dx dt = ax 1 − x K = ax − a K   x 2 (6)

Here, x is the size of the prey population at time t , a is an inherent per capita growth rate, and K is the carrying capacity. This model maintains a maximum prey population size, different from the original equations. When x = K, dx dt tends toward zero, and growth ceases, indicating that carrying capacity restricts indefinite growth.

For two prey populations, x 1 and x 2 , with logistic dynamics, the competitive Lotka-Volterra equations are: d x 1 dt = α 1 x 1 1 − x 1 + α 12 x 2 K 1 (7a) d x 2 dt = α 2 x 2 1 − x 2 + α 21 x 1 K 2 (7b) where α 12 represents the interactive effect species 2 has on the population of species 1 and α 21 represents the effect species 1 has on the population of species 2. All α -values are positive, indicating harmful interactions. Each species can have its own growth rate and carrying capacity (Bomze, 1995).

This model can be generalized for any number of competing species. Representing populations and growth rates as vectors and α 's as a matrix, the equation for any species i becomes (Equation 8): d x i dt = α i x i 1 − ∑ j = 1 N α ij x j K i (8)

When the carrying capacity is fully embedded within the interaction matrix, the resulting system (see Equation 9) modifies the standard expression accordingly: d x i dt = α i x i 1 − ∑ j = 1 N α ij x j (9) where N is the total number of interacting species. For simplicity, self-interacting terms α ii are often set to 1.

In adapting the Lotka-Volterra model to human demographic contexts, we consider regions as interacting entities whose populations influence one another through migration and other socio-economic factors (Wilson, 2008; Wilson and Dearden, 2011). We present a series of three equations representing two prey regions and one predator region, based on the foundations of the Lotka-Volterra model, with appropriate adjustments to account for the particular dynamics of human population evolution.

Each of the following nonlinear equations represents a specific geographical region.   X ˙ = AX − D X 2 + GXZ + G 1 XY  Y ˙ = CY − B Y 2 + G 4 YZ + G 2 XY  Z ˙ = − EZ + G 3 Z 2 + FXZ + G 5 YZ  (10)

These equations involve three categories of parameters derived and fine-tuned to reflect intrinsic growth or decline rates, as well as the interaction effects between regions X, Y, and Z.

The first category of computed parameters A, C, and E capture the intrinsic growth and decline rates of populations in regions X, Y, and Z, respectively. These parameters provide a baseline for intrinsic population dynamics within each region, representing how populations would behave under ideal conditions without external influences. They are computed using fixed parameters D, B, and G3, which are based on historical data and experiments. This ensures that the intrinsic growth and decline rates accurately reflect observed trends in each region.

Parameters A and C represent the intrinsic growth rates of prey regions X and Y, respectively, while parameter E reflects the intrinsic decline rate of predator region Z. They are derived using the formulas: A = Initial   Population X · D C = Initial   Population Y · B E = Initial   Population Z · G 3

The second category of D , B , and G 3 fixed parameters are crucial in population dynamics model, representing density-dependent effects within each region’s population. These parameters are fine-tuned experimentally to ensure the model accurately predicts populations for the years for which populations are already counted (known). The process involves iteratively adjusting the values of D, B, and G3 until the model’s predictions align very closely with actual population data. This fine-tuning going under an optimization process is essential for the robustness of the model in predicting known populations to which corresponding parameters’ values are derived, necessary for accurately predicting the future population. Parameter D represents the density-dependent factor for the population in region X, where higher values indicate stronger regulation, slowing growth as it nears carrying capacity. Similarly, parameter B affects region Y, with higher values implying more significant growth regulation. Parameter G 3 pertains to region Z, where increased values accelerate population decline as it approaches carrying capacity. These parameters account for limiting effects like resource and space constraints, ensuring the model reflects realistic population dynamics.

The second category includes the fixed parameters D, B, and G3, which represent density-dependent effects within each region’s population. These parameters are fine-tuned experimentally to ensure the model accurately predicts known populations. This process involves iteratively adjusting the values of D, B, and G3 until the model’s predictions closely align with actual population data. Such fine-tuning is essential for the robustness of the model in predicting future populations. Specifically, parameter D represents the density-dependent factor for the population in region X, where higher values indicate stronger regulation, slowing growth as it nears carrying capacity. Parameter BBB affects region Y in a similar manner, while parameter G3 pertains to region Z, where increased values accelerate population decline as it approaches carrying capacity.

The third category of interaction parameters G, G1, G4, G2, F, and G5 captures the dynamic interdependencies between the populations of regions X, Y, and Z. These parameters quantify how population changes in one region influence the others and are derived using the annual growth rate (AC) values of the respective regions, scaled by an adjusted gravity constant. By incorporating these AC values, the model accurately reflects the probabilistic and dynamic nature of population interactions. For instance, G and G1 capture the effects of region X on regions Z and Y, respectively, while G2 and G4 reflect the reciprocal influences between regions Y and X, and Y and Z. Similarly, F and G5 measure the impact of region Z on regions X and Y. Integrating these parameters into the adapted Lotka-Volterra equations ensures the model accurately simulates and predicts future population changes, providing a robust framework for understanding complex population dynamics and inter-regional influences.

The interaction parameters in the above triadic Equation 10 are calculated based on Newton’s formula. The gravity model in population dynamics is adapted from Newton’s law of gravity, where the force between two masses is proportional to their sizes and inversely proportional to the distance between them. In classical physics, Newton’s law of gravity states:

The interaction parameters in the above triadic Equation 10 are calculated based on Newton’s formula. The gravity model in population dynamics is adapted from Newton’s law of gravity, where the force between two masses is proportional to their sizes and inversely proportional to the square of the distance between them. In classical physics, Newton’s law of gravity is given by Equation 11, which states: F = K m 1 · m 2 d 2 (11) where F is the gravitational force, K is a fixed gravitational constant value, m 1 and m 2 are the masses of the two objects, and d 2 is the distance between them.

Recognizing the limitations of the traditional gravity model, Wilson (1970), Wilson (2010) introduced an entropy-maximizing approach to spatial interaction, providing a theoretical foundation that extends the gravity model by incorporating principles from statistical mechanics. Wilson’s gravity model for interregional flows can be written as Equation 12: T i j = O i D j   f c i j ∑ k D k f c i k (12) where:

• T i j is the flow from origin region i to destination region j .

Wilson’s model emphasizes that the interaction between regions depends not only on their masses (analogous to population sizes) and distances but also on the distribution of opportunities and the constraints of the system, ensuring that the flows are consistent with an entropy-maximizing principle.

In our approach, we adapt the gravity model to more accurately represent regional population interactions by ensuring that derived parameters closely reflect observed population dynamics. Unlike Wilson’s entropy-maximizing model, which determines the most probable distribution of flows under specific constraints, our model introduces three significant modifications to the classical gravity model. Additionally, we integrate the enhanced gravity model within the Lotka-Volterra framework, thereby improving its realism and applicability for population forecasting.

First, we replace the gravitational constant K with a probabilistic gravity weight Grv, reflecting the relative population sizes of the interacting regions. This change makes the model more dynamic and reflective of real-world population interactions, where the influence between regions is not static but varies with their population proportions. Grv i = P i P i + P j (13a) Grv j = P j P i + P j (13b)

Here, Grv i represents the probability of the population of region i relative to the total population of regions i and j . Similarly, Grv j represents the probability of the population of region j relative to the total population of regions i and j . Notably, Grv i ≠ Grv j , and Grv i + Grv j = 1 , contrary to the original fixed constant in Newton’s formula. For instance, in the differential Equation 10, the term GXZ represents the consumption of Region X by Region Z. The computation of the parameter G, in comparison to the parameter F in the term FXZ, will have different Grv effects. The calculation of Grv i ​ and Grv j aims to distinguish between consumption effects on the prey and nourishment effects on the predator, both subject to the probabilities of population.

Second, to simplify computations, we set d 2 to 1, assuming that distances between regions are uniform or negligible. However, it is possible to modify the distance term in various ways, ranging from simple to more complex specifications. For instance, Bernard et al. (2016) discussed the importance of refining distance measures in gravity models to better capture internal migration patterns, emphasizing how variations in distance decay can impact the accuracy of migration predictions. Similarly, Lewer and Van den Berg (2008) incorporated variations in socioeconomic distance, such as differences in income levels, to improve the explanatory power of gravity models in international immigration studies. These modifications illustrate the potential for fine-tuning the distance parameter to achieve more realistic modeling. Nonetheless, our study maintains d 2 = 1 to focus on the unified integration of the Lotka-Volterra and Gravity models, ensuring computational simplicity while effectively modeling population dynamics and spatial interactions.

Third, instead of using the product of the masses m 1 · m 2 , “population sizes” and their attractiveness or repulsiveness, as considered by authors like Curry (1972), Biagi et al. (2011), and Stillwell et al. (2014) or the attractiveness terms O i D j   in Wilson’s model, we utilize the annual growth rates AC i and AC j of populations in different regions. This approach provides a more accurate measure for several reasons:

• Annual growth rates account for temporal changes, reflecting current demographic trends rather than static population sizes. This aligns with the dynamic interactions in the Lotka-Volterra equations.

The annual growth rates are calculated as: AC i = ∆ P i P i (14a) AC j = ∆ P j P j (14b) where AC i represents the probability of growth of the population in region i ; AC j represents the probability of growth of the population in region j . ∆ P i is the population increase in region i and P i is the initial population in region i .

As illustrated by Equation 15, the interaction of two regions is the product of their AC and Grv probabilities. Excluding these factors reduces the model to traditional methods that rely solely on population size. The interaction parameters in the triadic Equation 10, such as G are thus derived using these probabilistic weights, providing a robust framework for demographic analysis and future population predictions: G = Grv i × AC i × 100 × AC j × 100 (15)

Here, G refers to the probabilistic rate of attraction (or repulsion) between two regions i and j , considering both the probabilistic constant Grv and the probabilistic growth rate of populations, both integrated within the same law of gravity. This Formula 14 exemplifies the calculation for other interaction parameters such as G1, G2, G3, G4, F, and G5.

To accurately estimate the parameters of the model and predict population dynamics, we employed an optimization process using Python programming. Specifically, we utilized the fsolve function from the scipy. optimize library to solve the dynamical system of equations. This function iteratively adjusts the values of X, Y, and Z to find the roots of the equations, effectively identifying the equilibrium points where the growth rates of the populations balance out. The optimization involves adjusting the parameters iteratively to minimize the error between the model’s predictions and observed data.

While acknowledging the substantial contributions of Wilson (1970), Wilson (2010) and others to spatial interaction models, this study diverges from traditional gravity-based frameworks by embedding them within a dynamic, probabilistic context. Earlier approaches often relied on static parameters and did not fully account for the nonlinear feedback mechanisms and evolving spatial interactions that characterize human population dynamics. By contrast, our method integrates a modified gravity model with Lotka-Volterra equations, incorporating probabilistic weights, dynamic growth rates, and carrying capacities. This unified framework captures the complexities of population mobility and interdependencies among regions, thereby offering a robust and forward-looking tool for demographic forecasting.

This innovative integration bridges the gap between ecological modeling and spatial interaction theories. Unlike previous theoretical studies that treated Lotka-Volterra or gravity models in isolation, we combine their strengths to produce more accurate and realistic predictions of population trends. The resulting model can simulate competitive and cooperative interactions among three regions and predict their movement toward an equilibrium state. Although extending the model to more regions is possible, it would introduce significant analytical complexity. For now, focusing on three regions provides a balanced, interpretable demonstration of the model’s capabilities.

By enabling not only population predictions but also the evaluation of equilibrium states, our approach increases both the robustness and comprehensiveness of demographic analysis. This improved representation of complex population behaviors supports more effective strategic urban and regional planning (Echenique, 1994). Going beyond static or deterministic perspectives, future developments could incorporate Bayesian inference to estimate parameters and yield probability distributions, further enhancing uncertainty quantification. In doing so, the model stands to advance the field of demographic forecasting and contribute valuable insights for policymakers and planners engaged in shaping sustainable demographic futures.

While this study’s primary innovation lies in adapting the Lotka-Volterra framework with a probabilistic gravity-based approach to better capture human population dynamics, future refinements could benefit from incorporating Bayesian inference as suggested by Wilson et al. (2021) and Raftery and Ševčíková (2023). Rather than altering the model’s foundational theoretical perspective, Bayesian methods would complement and refine its parameters, providing probability distributions and predictive intervals that more accurately reflect uncertainty and variability in demographic processes. By combining this strong theoretical core with advanced probabilistic inference, the model could achieve greater robustness and adaptability, offering a more nuanced and uncertainty-aware perspective on regional population forecasting.

In this section, we present a demonstration of population prediction based on the Lotka-Volterra differential equations. For this demonstration, we use data from the Insured Persons Registration File (IPRF) of the Health Insurance Authority of Quebec, collected and publicly shared by the Statistical Institute of Quebec. Our model utilizes population data incremented annually between 2021 and 2023 for regions X, Y, and Z, serving as the baseline for future population estimates. Additionally, we incorporate data on interregional migration, organized by age groups, to capture nuanced mobility patterns.

Our model is grounded in established theories from ecological modeling and spatial interaction, previously tested in various contexts (Wilson, 2008; Wilson and Dearden, 2011). We validated our model by comparing the predicted populations for 2021 to 2023 with actual historical data (Statistics Canada, 2023). The close alignment between predictions and observed values confirms the model’s accuracy and reliability. Detailed comparisons are provided in Table 2 and discussed below.

The analysis of Table 1 highlights specific regional gains and losses in 2021. For instance, Montreal experienced a population loss of 63,615, about 31% of all provincial outgoing movements. Of these, 41,633 people moved to its immediate periphery, represented by western regions including Laval, Lanaudière, and Montérégie. The central regions received 14,126 people, with Laurentides standing out, followed by Estrie and Mauricie. Another 4,836 moved eastward, concentrated around the Quebec City metropolitan area, with Outaouais receiving 1,340 people.

Conversely, Montreal gained 29,146 interregional immigrants in the same year, resulting in a net loss of 34,469 people. Specifically, 18,551 of these immigrants came from adjacent western cities. This figure, less than half of what Montreal lost, indicates the central area of Montreal is notably extending towards its peripheral crown. Immigration to the central region where Montreal is situated included 6,247 individuals from the central region and 2,939 from the eastern region.

Interprovincial migration trends since the 1960s indicate that Quebec has traditionally lost population to other provinces (Statistics Canada, 2010; Statistics Canada, 2013). However, over the past decade, this outflow has significantly slowed, nearly neutralizing the net balance. Data from the Ministry of Immigration and Cultural Communities show a consistent increase in net international migration to Quebec (Statistics Canada, 2023). Both interprovincial and international migration levels have remained relatively stable, often cancelling each other out. In contrast, interregional movements within Quebec strongly influence its population dynamics. Our forecasting model thus focuses on interregional mobility among the province’s 17 administrative regions.

As shown in Figure 1, to manage the complexity of multiple equations and enhance analytical clarity, we grouped the data into three main geographical divisions: the Montreal Metropolitan Region (RMR Montreal), the Quebec Metropolitan Region (RMR Quebec), and the Rest of the Province. These divisions are based on historical development, geographical proximity, and socio-economic similarities, providing a streamlined and insightful analysis. Among these three regions, RMR Montreal is considered the predator (Region West or Z), while the other two are considered prey (Region Center or X and Region East or Y). As indicated previously, unique conditions are integrated into each differential equation for these groupings.

We used Python software to codify parameters and conduct simulations to estimate the populations for regions X, Y, and Z for the years 2021–2023, as well as for projections from 2024 to 2053. Table 2 provides a summary of these results, including their evolutive parameter estimates, existing populations, and projected populations. The model predicts the population dynamics of regions X, Y, and Z over 30 years using the Lotka-Volterra equations and the estimated parameters based on the gravity model. The predicted populations are compared with known historical data to validate the model. The close alignment of predictions with known data indicates the model’s accuracy and robustness. This validation process ensures that the model reliably captures the population dynamics and interactions.

Also, with the increased population values in 2022 (reported on the populations of 2021 as initial values for 2022), parameters D, B, G3 and A, C, and E are recomputed. This allows us to predict known populations for 2021, 2022, and 2023 and derive all parameters and their trends over time. This iterative process of population and parameter forecasts continues yearly from 2024 to 2053.

The analysis of parameter trends provides a detailed understanding of population growth patterns across Quebec’s regions, indicating stabilization in growth rates and reduced intra-species competition. Table 2 presents the parameter values utilized in the model, along with their derivation and calibration, ensuring transparency and replicability. The results, supported by the data and graphics, highlight distinct population trends for regions X, Y, and Z from 2021 to 2053. For example, Region X is projected to grow from 1,611,613 in 2023 to 2,142,541 by 2053, with an average annual increase of approximately 17,698 individuals. This projection is consistent with historical growth trends for the region, which benefits from its strategic location between Montreal and Quebec City, attracting population influx due to its geographic and socio-economic advantages.

Region Y, encompassing the eastern part of the province, demonstrates a more modest increase in population. From 2023 to 2053, the population is projected to grow from 1,937,103 to 2,075,614 in Region Y, translating to an average annual increase of 4,617 individuals. Despite significant increases in previous years, the projection indicates stabilization in population growth, reflecting historical out-migration patterns where younger populations relocate to central and western regions.

Region Z, including the metropolitan area of Montreal, Laval, and Longueuil, shows an almost negligible increase in population from 5,224,710 in 2023 to 5,242,620 in 2053, with an average annual increase of 597 individuals. The decreasing growth trend for Region Z suggests that the population of the Montreal metropolitan area is gradually dispersing to peripheral regions, particularly to Region X, which offers better socio-environmental quality, less density and congestion, and more affordable living conditions (for instance, lower rents, housing prices, and taxes).

The analysis of population and parameter estimates provides a comprehensive understanding of their dynamic growth patterns in Quebec’s regions. Region X shows the most substantial growth, reflecting its strategic geographical location and socio-economic opportunities. Region Y demonstrates a stabilizing trend, indicating minimal growth due to out-migration. Region Z sees its population almost unchanged, potentially due to its population being drawn to peripheral regions. Graphic 1A corroborates these findings, showing convergence towards equilibrium points. This comprehensive analysis underscores the robustness of the Lotka-Volterra model in capturing complex population dynamics and providing reliable predictions for future population trends.

In comparison to our approach, the Cohort Component Model (CCM) results in Graphic 1B predict three different population growth scenarios for Quebec until 2056 (ISQ, 2019). The reference scenario projects the population to reach 9.9 million by 2056, with other scenarios estimating a range between 7.8 million (weak scenario) and 12 million (strong scenario). In comparison, our model, which incorporates Lotka-Volterra equations integrated with a probabilistic gravity model, predicts a total population of approximately 9.5 million by 2053. Both models show a growing trend, indicating that Quebec’s population is expected to increase in the coming decades. The similarity in trends between the two models demonstrates the robustness of our approach, which considers interactions and competition between regions, providing a more dynamic and nuanced prediction compared to the classic CCM approach. The additional insights from the ISQ (2013); ISQ (2019) study using the CCM approach highlight factors such as fertility, migration, and natural increase, which are also implicitly captured in our model through the interaction parameters and regional dynamics.

We have thoroughly defined and explained the parameters in the model specification section, including parameters such as A, C, E, D, B, and G3, along with interaction parameters like G, G1, and G2. These parameters form the foundation for the analysis of population growth trends in Quebec’s regions, and their detailed values are presented in Table 2. The analysis of these parameter trends suggests a shift from competitive to cooperative interactions over time, leading to stabilized growth rates and reduced intra-species competition. The results, as specified in Table 2 and further illustrated in the graphics, provide a comprehensive view of these population dynamics across several years, as projected by the Lotka-Volterra model.

The parameters D, B, and G3 show a general decline from 2021 to 2053. Specifically, D decreases from 1.732 in 2021 to 0.0080 in 2053, indicating a diminishing rate of intra-species competition or density-dependent effects within the populations. The parameter B, which relates to the intrinsic growth rate, also declines from 1.100 to 0.0700, suggesting a reduction in the natural growth rate over time. The parameter G3, representing inter-species interactions, transitions from a negative value (−0.552) to a slightly positive one (0.0300), implying a shift from competitive to more neutral or even cooperative interactions among the regions.

Parameters A, C, and E exhibit significant variability, with A and C decreasing and stabilizing around 0.3400 by 2053, while E shows a sharp reduction from −2.834 to 0.1400. This reflects changes in external factors such as immigration and economic influences affecting the growth dynamics of the regions.

The parameters G, G1, G4, G2, F, and G5 also show a decreasing trend, stabilizing towards the end of the period. This indicates that the inter-regional interaction coefficients are becoming less influential over time, possibly due to a more balanced distribution of population and resources among the regions.

The evolution of these parameters underscores the dynamic nature of population interactions and the gradual stabilization of regions. Over time, the parameters begin to stabilize, showing lesser variability, which implies that the model has reached a more robust and steady-state configuration. Parameters transitioning from positive to negative or approaching zero suggest a shift from strong interactions to more neutral or independent behaviors among regions. The fact that the parameters tend to linearize and stabilize suggests that the model’s predictive power is becoming more reliable as it adjusts to the underlying population dynamics. This comprehensive analysis provides valuable insights into the long-term behavior and stability of the population dynamics model, demonstrating how inter-regional interactions evolve over time and impact overall population predictions.

The three graphics collectively provide a comprehensive picture of the parameter estimates and population predictions over time for the XYZ regions. Graphic 2A illustrates the stabilization of parameter values over time, indicating reduced dynamic interactions and suggesting the system governed by the Lotka-Volterra model is reaching a quasi-equilibrium state. The lines representing the parameters converge and stabilize, reflecting the long-term balance achieved through inter-regional interactions and internal growth constraints.

Graphic 2B shows the population predictions for regions X, Y, and Z, with each region trending towards stabilization. Region Z, which initially shows significant fluctuations, eventually stabilizes. Regions X and Y also show trends towards stabilization, albeit at different rates. The total population curve exhibits a consistent growth trend, reflecting the stable and converging estimates over the long term. This stability can be attributed to the stabilization of parameter estimates, providing more consistent and reliable inputs into the Lotka-Volterra equations.

Graphic 2C presents a 3D view of the population predictions, showing how populations in regions X, Y, and Z evolve and converge over time. The spiraling behavior towards a central point indicates that the populations are interacting and adjusting to reach a stable equilibrium point by 2053. The result obtained in 3D clearly shows the convergence of X, Y, and Z from the wider periphery towards the center, reaching an equilibrium point, as evidenced by the spiral behavior observed. This convergence point is known as the equilibrium or prediction point, with values of 2.14 million for X, 2.07 million for Y, and 5.24 million for Z in 2053.

The convergence point estimate in 2053 underscores the robustness and validity of our model, effectively capturing the interdependencies and interactions between regions. It demonstrates that the populations of the three regions are not just stabilizing independently but are also interacting in a manner that leads to a collective equilibrium at the provincial scale. The model’s robust predictions highlight its ability to reflect complex population dynamics and accurately forecast long-term population trends. This interaction and convergence are driven by the probabilistic nature of the gravity model, which accounts for the relative attractiveness and repulsiveness of each region, resulting in a balanced and dynamic interaction.

By integrating advanced spatial interaction modeling techniques with the adapted Lotka-Volterra framework, our approach provides a robust and theoretically sound basis for predicting human population dynamics. It effectively addresses the complexities of migration patterns and spatial interactions through an improved probabilistic gravity model. The predictions from this unified approach closely align with observed data, confirming both its reliability and applicability to demographic forecasting. For instance, the predicted population for Region X in 2023 was within 1% of the actual recorded population, with similarly high accuracy observed for Regions Y and Z. This calibration not only validates the model’s theoretical assumptions but also ensures its practical utility for long-term population projections by accurately reflecting observed trends.

The convergence point itself is a conclusive and highly satisfactory result, as it could have diverged towards non-equilibrium points, yielding abstract and unrealistic outcomes. The equilibrium shows a “demographic attractive point” towards which the system evolves, similar to other complex systems like the Lorenz (1963) attractor arising in a system of dynamic equations under other types of forces (temperature difference, gravity, buoyancy, kinematic viscosity, etc.). That point can be located in space as an epicenter in Quebec province in 2053, allowing planners and policymakers to take strategic decisions on resource planning and sharing.

Accurate demographic forecasting is essential for strategic urban planning and effective resource allocation, particularly in regions characterized by complex, nonlinear population dynamics shaped by socioeconomic factors and spatial interactions. Traditional methods, such as the Cohort Component Model (CCM), tend to rely on linear projections, which often fail to capture the intricate interdependencies driving population fluctuations.

To address these limitations, this study presents an innovative approach that integrates the Lotka-Volterra prey-predator model—originally designed for ecological interactions—and adapts it for human demographic forecasting, thus providing the theoretical and conceptual foundations often lacking in traditional analyses. This adaptation is further strengthened by incorporating a rigorously specified probabilistic gravity model, uniquely embedded within the Lotka-Volterra framework. By merging these approaches, we account for both competitive and cooperative dynamics shaping population growth, as well as the spatial mobility between regions.

The application of the Lotka-Volterra model to human demographic contexts enables a detailed representation of nonlinear population dynamics and predictions. By incorporating a modified gravity model—rooted in spatial interaction theories—this study considers migration flows governed by geographic proximity and socioeconomic attractiveness, offering a more refined approach to understanding population mobility. This integrated framework extends the existing literature in spatial interaction and urban systems modeling, providing a more comprehensive method for predicting regional population changes.

The model’s robustness and predictive accuracy are supported by empirical findings. The close alignment between predicted and observed population trends from 2021 to 2023, including a prediction error of less than 1% for Region X in 2023, attests to its reliability. Similarly high accuracy levels were observed for Regions Y and Z. These results confirm that the model effectively captures the complex dynamics driving population changes in Quebec. The observed shift from competitive to cooperative interactions in the parameter trends suggests stabilizing growth rates and regional equilibrium—a hallmark of the Lotka-Volterra framework—highlighting its applicability to demographic forecasting.

Although this study focuses on three regional clusters in Quebec, its theoretical and practical foundations are transferable to other geographic contexts. By extending the Lotka-Volterra model beyond its ecological roots and integrating it with the gravity model, this flexible framework offers a valuable tool for analyzing regional population dynamics in diverse settings. The methodology provides actionable insights for strategic planning, addressing the shortcomings of traditional models and delivering more accurate, long-term projections of population changes. By combining the Lotka-Volterra and gravity models, this approach yields more realistic predictions of population trends and spatial interactions.

The findings also highlight the model’s ability to identify distinct population growth trends across regions: Region X is projected to experience significant growth, driven by its strategic location and socioeconomic opportunities; Region Y shows more modest growth, influenced by historical out-migration; and Region Z, encompassing the Montreal metropolitan area, exhibits minimal growth, reflecting population dispersal to peripheral areas. These trends underscore the model’s ability to capture interregional interactions and offer critical insights for policymakers and urban planners, who can use these predictions to inform infrastructure development, resource allocation, and regional development strategies.

As with any complex model, there are opportunities for further refinement. Incorporating additional variables, such as economic indicators, policy changes, and more nuanced measures of distance, could further enhance predictive accuracy. Moreover, integrating micro-level dynamics, such as individual decision-making and migration motivations, may improve the granularity of future forecasts. Exploring advanced computational techniques, like agent-based modeling or machine learning, could also increase the model’s adaptability and precision, responding effectively to shifting demographic and socioeconomic conditions. These advancements contribute to more sustainable urban and regional development strategies, underscoring the importance of innovative approaches in meeting the challenges of modern demographic forecasting.

Looking ahead, while the current integration of Lotka-Volterra dynamics with a probabilistic gravity model provides a theoretically grounded, equilibrium-oriented forecasting framework, it also suggests new avenues for refinement. At present, the model’s predictive capacity relies on systematically adjusted parameters and conditions guiding it toward realistic equilibria. Future research could incorporate explicit probabilistic methods—such as Bayesian inference—to better represent uncertainties in parameters and evolving socioeconomic landscapes. Additionally, exploring complementary methodologies, including entropy-based measures, could enrich analytical depth. Rather than replacing the model’s core theoretical insights, these enhancements would refine its ability to accommodate variability and change over time, ultimately offering a richer understanding of complex demographic processes and enabling more adaptable, regional planning–relevant forecasting.