Mathematics and Computers in Simulation 234 (2025) 113–134

Contents lists available at ScienceDirect

Mathematics and Computers in Simulation

journal homepage: [Link]/locate/matcom

Original articles

An agile modeling framework for population dynamics

Laurent Attias a ,∗, Vincent Siess b , Stéphane Labbé a
a
Sorbonne Université, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions, Paris, F-75005, France
b
Dassault Systèmes, 10 rue Marcel Dassault, Vélizy-Villacoublay, 78140, France

ARTICLE INFO ABSTRACT

MSC: In this work, we present an agile modeling framework for structured population dynamics,
92D25 leading to automated generation of population models’ equations. The structure of a population,
92-10 i.e. its separation in strata, according to one or many criteria (such as sex, income, health,
35Q49
geographic area or species if dealing with animal populations), represents a major issue for the
34K20
precision and richness of population dynamics simulations. The intensity of some phenomena
92D30
65P30
and mechanisms is highly dependent on the involved subpopulation characteristics. This mod-
eling framework can be seen as an extension of the classic McKendrick–von Foerster equation,
Keywords:
which embeds the population structure. It allows showing, under appropriate hypothesis, an
Population dynamics
McKendrick equation existence and local uniqueness result for the solution of a transport equation. A modeler
Structured population models has been implemented, to generate models that respect the desired structure hypotheses. We
Predator–prey model illustrate its abilities on an age-structured predator–prey model, subject to migratory dynamics
Model generation and to an epidemic, based on a SIRD model.
Bifurcations

0. Introduction

Modeling and simulating population is a major issue in our modern societies. Forecasting population movements, changes in
composition, health evolutions represent key features for governments, public institutions or private insurances. In this work, we aim
to provide a unified framework for population dynamics modeling and simulation. As age is an important feature, we generalize the
classic McKendrick–von Foerster equation, in order to take into account various phenomena in a multi-structured population, such
as exchanges between compartments. These are defined as subpopulations sharing common traits. We will rely on this framework to
present a modeler, i.e. a software able to generate simulable models in a high-level model description language like Modelica. Our
goal is to be able to model and simulate any structured population problem. For example, studying the dynamics of a structured
Lotka–Volterra model, including aging, migration and an epidemic, would be a satisfying objective.
Since the end of the 18th century, modeling a population with discrete or differential equations has been well known and has
proved its efficiency in a first approach [1]. But to get a precise simulation, heterogeneity needs to be taken into account, which
leads to structured models of population. Such models describe a population as a set of subpopulations, defined with one or several
traits (sex or gender, size, age, geographic area...). These traits which may influence natality, growth or mortality rates, as well
as interactions between individuals within a compartment or between compartments. The study of structured population aims to
analyze how those traits affect the model dynamics [2]. Structured population models are highly used in biological mathematics,
including the study of cells growth [3] or epidemic dynamics [4]. To model heterogeneity and provide insights on differences
between individuals, stage-structured models have been developed, sometimes within a non-deterministic point of view [5]. Related

∗ Corresponding author.
E-mail addresses: [Link]@[Link] (L. Attias), [Link]@[Link] (V. Siess), [Link]@[Link] (S. Labbé).

[Link]
Received 22 July 2024; Received in revised form 13 December 2024; Accepted 10 February 2025
Available online 21 February 2025
0378-4754/© 2025 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights are
reserved, including those for text and data mining, AI training, and similar technologies.
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

works deal with life expectancy through stage and age-structured models [6]. Structured populations may also be studied from the
genetic point of view [7,8], but this is not the approach retained in this work.
The McKendrick–von Foerster equation [9] is one of the most fundamental equations of population dynamics. This transport
equation describes the evolution in time of an age-structured population, with births, aging and death. This model appeared for the
first time in the context of epidemiology [10,11], but its thorough analysis was achieved in a later step [12,13], thanks to methods
related to Volterra integral equations and Laplace transformation.
More recently, the introduction of non-linear models breathed new life into age-structured populations study. In Gurtin and
Maccamy [14], using non-linear integral Volterra equations, authors showed the existence, uniqueness and convergence to an
equilibrium of solutions to a non-linear Sharpe–Lotka–McKendrick model. This work paved the way to breakthroughs on non-linear
models, as much on theoretical developments as on biological applications (see [15] for an exhaustive approach). The growing
mathematical complexity of non-linearities in age-structured models led to developments of new tools and methods, such as linear
and non-linear operator semigroups in Banach spaces, with a functional analysis approach [16].
Structured population models has been studied in a cumulative formulation [17], which consists in considering various
mechanisms (aging, reproduction, death, migration. . . ) separately and adding their contributions, leveraging the superposition
principle on linear differential equations. In a such approach, population dynamics can be studied at both individual and population
level (𝑖-level and 𝑝-level, respectively) [18]. At the 𝑖-level, a probabilistic point of view prevails, as equations such as the Chapman–
Kolmogorov involve the probability to reach a state in a given subset of the states space, for an individual at a given state, at
a given time. At the 𝑝-level, a book-keeping operates to aggregate changes at the 𝑖-level. The renewal of the population, due to
reproduction and deaths, is expressed through integral equations that can be related, if one considers age as a transport variable,
to McKendrick-like equations.
The branching processes field introduced the terminology of kernel [19,20] to describe an elementary operator which rules a
phenomenon’s dynamic (e.g. the reproduction kernel).
In their works, O. Diekmann, M. Gyllenberg, J.A.J. Metz and H.R. Thieme built a strong formal and theoretical framework,
enabling to derive mathematical properties from population models (see Diekmann et al. [17,18,21,22]). In this paper, our goal is
different and consists in building a formal framework to enable a consistent numerical implementation.
The McKendrick–von Foerster equation is still used in both epidemiology and demography studies [23]. For further reading on
transport equations and application to structured population modeling, one should consult [24].
At Dassault Systèmes, we have been interested in population dynamics simulation for many years. We used an algebraic–
differential equation (DAE) formalism, representing subpopulations as compartments, relying on the Modelica language [25], a
high-level model description language. Solvers such as Dymola [26] lie on the Modelica language to perform numerical integration
of DAEs. Our objective is to develop a population model and simulation that could be adapted and tuned to stakeholder’s hypotheses
and needs, for instance in terms of structuration, granularity, scale, or scope of the simulation. Unfortunately, editing manually huge
Modelica models, in order to change its structuration and behavior laws, is time-consuming, painful and presents an important
error risk. There was a need for a tool able to take as input comprehensive population structuration and demographic phenomena
description, and to provide as output an executable file which will then be given to the solver for simulation.
To support such generality, this tool needed a unified framework – a metamodel – able to take into account any population
structure and describe any demographic phenomenon. The McKendrick–von Foerster equation turned out to be a remarkable
complete metamodel for this objective. It can be instantiated with a DAE model, that is what we did, but also within other paradigms,
such as Monte-Carlo or multi-agent simulations.
This article is divided in three sections. In the first section, we present a modeling framework for structured populations, based
on an extension of the McKendrick–von Foerster equation. In the second section, we present a model generator developed and
implemented to generate structured population models based on the previous formalism. In the last section, we apply this formalism
and this model generator to a structured Lotka–Volterra model, and briefly discuss numerical results.

1. Modeling framework

1.1. Extension of McKendrick formalism

1.1.1. The McKendrick–von Foerster equation

The McKendrick–von Foerster equation is a transport equation, highly used in biological and ecological mathematics. It models
the evolution in time of an age-structured population undergoing natality and mortality. The specificity of this equation compared to
a usual transport equation lies in its boundary condition: the births depend on the whole population. The original McKendrick–von
Foerster equation [10] is stated as follows:
⎧ 𝜕𝑡 𝜙(𝑡, 𝑎) + 𝜕𝑎 𝜙(𝑡, 𝑎) + 𝜇(𝑎)𝜙(𝑡, 𝑎) = 0 ∀(𝑡, 𝑎) ∈ R∗+ × R∗+ , (a)
⎪ +∞
⎨ 𝜙(𝑡, 0) = ∫ 𝑏(𝑎)𝜙(𝑡, 𝑎) d𝑎 ∀ 𝑡 ∈ R∗+ , (b) (1)
⎪ 0
⎩ 𝜙(0, 𝑎) = 𝜙0 (𝑎) ∀ 𝑎 ∈ R∗+ , (c)
where 𝜙 is the population density, depending on the time 𝑡 and the age 𝑎, 𝜇 is the mortality rate, depending on 𝑎, 𝑏 the birth rate,
also depending on 𝑎. It is generally assumed that 𝑏 has compact support in R∗+ , which ensures that 𝜙(𝑡, 0) can be computed from
[ ]
𝜙(𝑡, 𝑎) for 𝑎 in 𝑎min , 𝑎max with 0 < 𝑎min < 𝑎max . (1a)–(1c) can be solved in a closed form thanks to the method of characteristics
(see Murray [9]).

114
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

1.1.2. Generalization in a multi-structured context

We aim to generalize this equation in a unified formalism to deal with structured populations, with any structure one can imagine
(any number of traits, any number of values of these traits, continuous or discrete traits...). In this way, we consider the following
generalized McKendrick equation:
⎧ 𝜕𝑡 𝜙(𝑡, 𝑎, 𝑥) + 𝜕𝑎 𝜙(𝑡, 𝑎, 𝑥) + 𝐾(𝑡, 𝜙(𝑡), 𝑎, 𝑥) = 0 ∀ (𝑡, 𝑎) ∈ 𝐼 × R∗+ , for a.e. 𝑥 ∈ D, (a)
⎪ +∞
⎨ 𝜙(𝑡, 0, 𝑥) = ∫ 𝑏(𝑡, 𝜙(𝑡), 𝑎, 𝑥, 𝑦)𝜙(𝑡, 𝑎, 𝑦) d𝑦 d𝑎 ∀ 𝑡 ∈ 𝐼 , for a.e. 𝑥 ∈ D, (b) (2)
⎪ 0 ∫D
⎩ 𝜙 (0, 𝑎, 𝑥) = 𝜙0 (𝑎, 𝑥) ∀𝑎 ∈ R∗+ , for a.e. 𝑥 ∈ D, (c)
where:

∏
𝑁
• 𝑥 belongs to a set D called the hyperparameters set; typically, one could have D = D𝑖 where 𝑁 is the number of traits1 and
𝑖=1
D𝑖 the set of values of the 𝑖th trait.
• 𝐼 is an interval of R, of the form [0, 𝑇 ] or [0, +∞[,
( )
• 𝐾 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D → R is a global operator, defined for a.e. (𝑎, 𝑥) in R+ × D,
( )
• 𝑏 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D × D → R+ is the natality rate function, also defined for a.e. (𝑎, 𝑥) in R∗+ × D. 𝑏(𝑡, 𝑔 , 𝑎, 𝑥, 𝑦) represents
the natality rate at time 𝑡, for a population distribution at instant 𝑡 defined by 𝑔, for parents with age 𝑎 and hyperparameter
𝑦 giving birth to children with hyperparameter 𝑥,
( )
• 𝜙0 is the initial population, belonging to 𝐿1 R+ × D ,
( )
• 𝜙 is the population, considered at least in 𝐿1 𝐼 × R+ × D . We will later give a precise meaning of the partial derivatives in
(2a).
In Section 3.2.1, we present a concrete example of how this formalism can be applied to the modeling of a predator–prey situation
with age, space and epidemic structuration. For now, let us precise that the set of hyperparameters D is meant to be the product
of multiple structuration dimensions, representing continuous or discrete criteria. For example, we consider a population model
with, in addition to age, space, gender and income structuration. The age structuration is handled by the transport trait, so D will
gather only the other traits. If we consider 𝑛 distinct zones of space, two genders {𝖬, 𝖥}, and annual income in R+ , we will have
D = [[1, 𝑛]] × {𝖬, 𝖥} × R+ . A typical hyperparameter would be 𝑥 = (𝑛, 𝖥, 𝑟) to denote women in zone 𝑛 with annual income 𝑟. The
income structure can be described continuously, as in the previous example, but also discretely, replacing R+ by income groups
{ } [ [
𝐼1 , … , 𝐼𝑝 , where 𝐼𝑘 = 𝑖𝑘 , 𝑖𝑘+1 for 𝑘 in [[1, 𝑝]], with 𝑖1 = 0 and 𝑖𝑝+1 = +∞. The integral over D in (2b) is then converted into a
discrete sum on 𝑘 going from 1 to 𝑝, which is equivalent to considering the counting measure on the income criterion for integration
over D.
In (2b) we assume that the natality rate for a parent with (𝑎, 𝑦) characteristics giving birth to a newborn with 𝑥 characteristics
may depend on the whole population 𝜙(𝑡). For example, taking into account resources or influence between different compartments.
This is a first step of complexity compared to the original McKendrick boundary condition (1b). But we also assume that this natality
rate only depends on the population 𝜙(𝑡) at instant 𝑡, and not on the population evolution history (𝜙(𝑠))𝑠∈[0,𝑡] . This hypothesis is done
for the sake of simplicity and does not seem to be limitative from a modeling point of view, since the natality rate is at first order
influenced by real-time factors, such as economic development, access to education, status of women in society, access to healthcare,
and so on. This consideration also applies to global operator 𝐾, which may include mortality, geographic migration, evolution of
the income level, etc.
The combination of structuration criteria is simplified by the genericity of the presented framework. One could easily encapsulate
in the 𝑥 hyperparameter, belonging to the hyperparameter space D, each structuration criterion, such as geographic location, income
level or health status. The flows and relationships involved by the addition of a new criterion are easily taken into account in this
unified formalism.

1.1.3. Kernel and pointwise operators

As a global operator 𝐾, we may consider two subtypes of operators that will be useful for population dynamic modeling: pointwise
and kernel operators. The term of kernel can be linked to the terminology of Jagers [19,20] and the field of branching processes.

Pointwise operators. A pointwise operator describes an incoming or outcoming flow proportional to a subpopulation, whose
intensity may vary depending on the characteristics of the subpopulation but also on the global population at time 𝑡.
( )
We define it as a mapping 𝐾 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D → R, defined only for a.e. (𝑎, 𝑥) in R+ × D, so that for any (𝑡, 𝑔) in
1
( )
𝐼 × 𝐿 R+ × D and for a.e. (𝑎, 𝑥) in R+ × D,
𝐾(𝑡, 𝑔 , 𝑎, 𝑥) = 𝜇(𝑡, 𝑔 , 𝑎, 𝑥)𝑔(𝑎, 𝑥),
( )
where 𝜇 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D → R is a rate on R+ × D i.e. a mapping defined only for almost every (𝑎, 𝑥) in R+ × D; 𝜇(𝑡, 𝑔 , 𝑎, 𝑥)
represents the intensity of an incoming or outcoming flow within the (𝑎, 𝑥) compartment.

1 Excepted age, which is considered separately because it is the transport trait.

115
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Kernel operators. A kernel operator quantifies precisely the transfer flow to a given compartment from other ones. It is defined as
( ) ( )
a mapping 𝐾 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D → R, defined only for a.e. (𝑎, 𝑥) in R+ × D, so that for any (𝑡, 𝑔) in 𝐼 × 𝐿1 R+ × D and for
a.e. (𝑎, 𝑥) in R+ × D,
𝐾(𝑡, 𝑔 , 𝑎, 𝑥) = 𝜅 (𝑡, 𝑔 , 𝑎, 𝑥, 𝑦) 𝑔(𝑎, 𝑦) d𝑦,
∫D
( )
where 𝜅 is a kernel on R+ × D, i.e. a mapping 𝜅 ∶ 𝐼 × 𝐿1 R+ × D × R+ × D × D → R defined only for almost every (𝑎, 𝑥, 𝑦) in
R+ × D × D; 𝜅(𝑡, 𝑔 , 𝑎, 𝑥, 𝑦) represents the intensity of an exchange flow2 from the (𝑎, 𝑦) to the (𝑎, 𝑥) state in a population distribution
𝑔, at time 𝑡. It is assumed that a kernel may model a change of hyperparameter 𝑦 to 𝑥, but does not affect age, which evolves only
through the transport term.

1.1.4. Scope of the extended McKendrick formalism

From the two elementary types of operators, it is interesting noting that one can build conservative operators representing
( )
exchanges between compartments. Let 𝑘 be a rate of transfer, defined on 𝐼 × 𝐿1 R+ × D × R+ × D × D, defined only for almost every
(𝑎, 𝑥, 𝑦) ∈ R+ × D × D. We build the following operator:
𝐾(𝑡, 𝑔 , 𝑎, 𝑥) = (𝑘(𝑡, 𝑔 , 𝑎, 𝑦, 𝑥)𝑔(𝑎, 𝑥) − 𝑘(𝑡, 𝑔 , 𝑎, 𝑥, 𝑦)𝑔(𝑎, 𝑦)) d𝑦
∫D

We may see this exchange operator as a combination of a pointwise operator with rate (𝑡, 𝑔 , 𝑎, 𝑥) ↦ 𝑘(𝑡, 𝑔 , 𝑎, 𝑦, 𝑥) d𝑦 and a kernel
∫D
operator with kernel 𝑘. As we have
𝐾(𝑡, 𝑔 , 𝑎, 𝑥) d𝑥 = 0,
∫D
we recognize a conservative property of such operator; for example, conservation of the total population considering migrations
between territories.
The population modeling formalism presented is able to handle pointwise, kernel and conservative exchange operators, which
may be sufficient for a large range of applications. We acknowledge that for now, it cannot handle some types of operators which
lead to wider types of differential equations. The main limitation is that the global operator 𝐾 must depend on the population 𝜙
only through the current population 𝜙(𝑡). This point is a key feature to enable what will follow. Thus, delay equations involving
𝜙(𝑡 − 𝜏) or equations requiring a ‘‘memory’’ on 𝜙 with a dependence of 𝐾 on 𝜙(𝐴) for a subset 𝐴 of [0, 𝑇 ] are neither handled in the
extended McKendrick formalism.
A second limitation lies in modeling discontinuous features such as hybrid models involving discrete events, because of regularity
assumptions detailed in the following.
Some of these limitations, such as delay equations, can be circumvented thanks to transport equation formulations, that would
be compatible with both the modeling framework and the simulation method presented.
There exists a variety of models that can be embedded in the extended McKendrick formalism. Being able to build conservative
exchange operators paves the way to model a wide range of behaviors. For example, as we will develop below, the predator–prey
model and the Lotka–Volterra equations can be rewritten in this framework, same as epidemic compartmental models (SIR, SIRD,
and all their variations...). Modeling economic evolutions in a population, or information propagation, would also be possible within
our formalism.

1.2. Existence and local uniqueness

1.2.1. Fixed point formulation

We firstly reformulate (2a)–(2c) with 𝐿1 (D) objects, to get closer to (1a)–(1c) formulation:
⎧ 𝜕𝑡 𝜙(𝑡, 𝑎) + 𝜕𝑎 𝜙(𝑡, 𝑎) + 𝐾(𝑡, 𝜙(𝑡), 𝑎) = 0 ∀ (𝑡, 𝑎) ∈ 𝐼 × R∗+ , (a)
⎪ +∞
⎨ 𝜙(𝑡, 0) = ∫ 𝑏(𝑡, 𝜙(𝑡), 𝑎, ⋅, 𝑦)𝜙(𝑡, 𝑎, 𝑦) d𝑦 d𝑎 ∀ 𝑡 ∈ 𝐼, (b) (3)
⎪ 0 ∫D
⎩ 𝜙 (0, 𝑎) = 𝜙0 (𝑎) ∀ 𝑎 ∈ R∗+ . (c)

To give a precise meaning to the 𝜕𝑡 + 𝜕𝑎 transport operator, even if 𝜙 is not smooth, we reformulate (3a)–(3c) into a fixed point
formulation which will not involve such partial derivative operator, and will be equivalent to (3a)–(3c) if 𝜙 is smooth enough. This
reformulation is based on the characteristics method [9]. We define a mapping 𝑓 by
𝑓 ∶ (𝑡, 𝑎) ↦ 𝜙(𝑡, 𝑎 + 𝑡).

for any real 𝑎, not necessarily positive, and 𝑡 in the time interval [0, 𝑇 ]. This 𝑎 will be called pseudo-age, as 𝑎 + 𝑡 is the real age
considered and must be positive. We denote by 𝜙(𝑡) the mapping 𝑎 ↦ 𝜙(𝑡, 𝑎) for any age 𝑎. We then have 𝜙(𝑡) = 𝑓 (𝑡, ⋅ − 𝑡) and (3a)
rewrites, for fixed real pseudo-age 𝑎,
d𝑡 𝑓 (𝑡, 𝑎) = −𝐾 (𝑡, 𝑓 (𝑡, ⋅ − 𝑡), 𝑎 + 𝑡).

2 observe the convention on the order of origin 𝑦 and destination 𝑥 hyperparameters.

116
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Integrating this equation, distinguishing cases 𝑎 > 0 and −𝑡 ⩽ 𝑎 ⩽ 0 and taking into account initial and boundary conditions leads
to the following integral formulation, for all 𝑡 ∈ [0, 𝑇 ] and for a.e. 𝑎 ∈ [−𝑡, +∞[,
| 𝑡
| 𝜙 (𝑎) − 𝐾 (𝑠, 𝑓 (𝑠, ⋅ − 𝑠), 𝑎 + 𝑠) d𝑠 if 𝑎 > 0,
| 0 ∫0
| (
| +∞
|
| 𝑏 (|𝑎| , 𝑓 (|𝑎| , ⋅ − |𝑎|), 𝜏 , ⋅, 𝑦) 𝜙0 (𝜏 − |𝑎| , 𝑦)
| ∫|𝑎| ∫D
| )
𝑓 (𝑡, 𝑎) = | |𝑎| (4)
|
| − 𝐾 (𝑠, 𝑓 (𝑠, ⋅ − 𝑠), 𝜏 − |𝑎| + 𝑠, 𝑦) d𝑠 d𝑦 d𝜏
| ∫0
|
|
| 𝑡
| − 𝐾 (𝑠, 𝑓 (𝑠, ⋅ − 𝑠), 𝑎 + 𝑠) d𝑠 if −𝑡 ⩽ 𝑎 ⩽ 0.
| ∫
| |𝑎|

This formulation requires less regularity than (3a)–(3c) to be studied, and then may allow less regular solutions. One could think
that such integral formulation is making a mountain out of a molehill, as we could study regular solutions of (3a)–(3c), but there
is nothing of the sort. It is well known that in a transport problem, the regularity of the initial condition is usually reported to
the global solution. In population dynamics, we may imagine violent phenomena such as demographic shocks, wars, pandemics,
baby-booms, etc. Such events may be modeled as non-regular initial conditions. This justifies the study of a non-regular formulation
of the extended McKendrick equation.

1.2.2. Functional framework

To study (4), we introduce, for a time limit 𝑇 > 0, the following time-pseudo-age space:
𝑈𝑇 = {(𝑡, 𝑎) ∈ [0, 𝑇 ] × R, 𝑎 + 𝑡 ⩾ 0} .

We will look for solutions of (4) that will be measurable on 𝑈𝑇 , defined everywhere in time,3 𝐿1 (D) valued, for which the norm
defined by:
+∞
‖𝑓 ‖∞,𝐿1 = sup ‖𝑓 (𝑡, 𝑎)‖𝐿1 d𝑎
𝑡∈[0,𝑇 ] ∫−𝑡
( )
is finite. We will denote 𝐿∞,1 𝑈𝑇 the Banach space of such functions 𝑓 with ‖𝑓 ‖∞,𝐿1 < +∞.
We also define a specific norm on D × D:

Definition 1 (‖⋅‖𝐿1 -𝐿∞ Norm). Let 𝑢 be a function defined almost everywhere on D × D. We will denote by ‖𝑢‖𝐿1 -𝐿∞ the (potentially
infinite) quantity:
‖𝑢‖𝐿1 -𝐿∞ = sup ess |𝑢(𝑥, 𝑦)| d𝑥 = ‖ ‖
‖𝑦 ↦ ‖𝑢(⋅, 𝑦)‖𝐿1 ‖𝐿∞ .
𝑦∈D ∫D
( )
We will denote by 𝐿1 -𝐿∞ (D × D) the space of functions 𝑢 defined almost everywhere on D × D such as ‖𝑢‖𝐿1 -𝐿∞ < +∞.
( ) ( ) ( )
For 𝑟 > 0 and 𝑓0 in 𝐿∞,1 𝑈𝑇 , we will denote by 𝐵 𝑓0 , 𝑟 the closed ball with center 𝑓0 and radius 𝑟 in the space 𝐿∞,1 𝑈𝑇
endowed with ‖⋅‖∞,𝐿1 norm.

1.2.3. Assumptions on 𝐾 and 𝑏

We will assume the following properties on the global operator 𝐾 and the natality rate 𝑏.

Assumptions on the global operator 𝐾. We assume that:

(a) 𝐾 is integrable on R∗+ , in the sense that

( ) ( )
∀ (𝑡, 𝑔) ∈ 𝐼 × 𝐿1 R+ × D , 𝐾(𝑡, 𝑔) ∈ 𝐿1 R+ × D . (5)

( ) ( )
(b) 𝐾 is locally semi-Lipschitz, in the sense that for any 𝑡0 in 𝐼 and 𝑔0 in 𝐿1 R+ × D there exists 𝑉0 a neighborhood of 𝑡0 , 𝑔0
( ) ( ) ( )
inside 𝐼 × 𝐿1 R+ × D , and 𝑐0 ⩾ 0, such that for any 𝑡, 𝑔1 and 𝑡, 𝑔2 in 𝑉0 ,
‖ ( ) ( )‖
‖𝐾 𝑡, 𝑔1 − 𝐾 𝑡, 𝑔2 ‖ 1 ⩽ 𝑐0 ‖ ‖
‖𝑔1 − 𝑔2 ‖𝐿1 . (6)
‖ ‖𝐿

( )
(c) 𝐾 is bounded on every compact, in the sense that for any compact subset 𝑉 of 𝐿1 R+ × D , there exists 𝑀𝐾 ,𝑉 ⩾ 0, such
that for any (𝑡, 𝑔) in 𝑉 ,
‖𝐾(𝑡, 𝑔)‖𝐿1 ⩽ 𝑀𝐾 ,𝑉 . (7)

3 a function 𝑓 defined on 𝑈 will be said defined everywhere in time if it is defined almost everywhere on 𝑈 , with for any 𝑡 in [0, 𝑇 ], 𝑓 (𝑡, ⋅) ∶ [−𝑡, +∞[→ 𝐿1 (D)
𝑇 𝑇
is well-defined.

117
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Assumptions on the natality rate 𝑏. We assume that:

(a) 𝑏 is 𝐿1 -𝐿∞ -integrable, in the sense that

( ) ( )
∀ (𝑡, 𝑔 , 𝑎) ∈ 𝐼 × 𝐿1 R+ × D × R+ , 𝑏 (𝑡, 𝑔 , 𝑎) ∈ 𝐿1 -𝐿∞ (D × D) (8)

( uniform
(b) 𝑏 has ) compact support in age 𝑎. More precisely,
[ there
] exists 0 < 𝑎min < 𝑎max such that for any 𝑡 in 𝐼 and 𝑓 in
𝐿1 R+ , 𝐿1 (D) , the support of 𝑏(𝑡, 𝑓 ) is included in 𝑎min , 𝑎max .
( )
(c) 𝑏 is essentially uniformly 𝐿1 -𝐿∞ -bounded, in the sense that for any 𝑡 in 𝐼 and 𝑔 in 𝐿1 R+ × D ,
( ( ) )
𝑏(𝑡, 𝑔) ∈ 𝐿∞ R+ , 𝐿1 -𝐿∞ (D × D) , (9)
where ‖𝑏(𝑡, 𝑔)‖𝐿∞ = sup ess ‖𝑏(𝑡, 𝑔 , 𝑎)‖𝐿1 -𝐿∞ .
𝑎∈R
+
1
( ) ( )
(d) 𝑏 is locally semi-Lipschitz,
( ) in the sense that for any 𝑡0 in 𝐼 and
( 𝑔)0 (in 𝐿 ) R+ × D , there exists 𝑉1 a neighborhood of 𝑡0 , 𝑔0
1
inside 𝐼 × 𝐿 R+ × D , there exists 𝑐1 ⩾ 0, such that for any 𝑡, 𝑔1 , 𝑡, 𝑔2 in 𝑉1 ,
‖ ( ) ( )‖
‖𝑏 𝑡, 𝑔1 − 𝑏 𝑡, 𝑔2 ‖ ∞ ⩽ 𝑐1 ‖ ‖
‖𝑔1 − 𝑔2 ‖𝐿1 . (10)
‖ ‖𝐿

(e) 𝑏 is uniformly 𝐿1 -𝐿∞ -bounded, in the sense that there exists 𝐵 > 0 such that
( )
∀ (𝑡, 𝑔) ∈ 𝐼 × 𝐿1 R+ × D , for a.e. 𝑎 ∈ R+ , ‖𝑏 (𝑡, 𝑔 , 𝑎)‖𝐿1 -𝐿∞ ⩽ 𝐵 . (11)

1.2.4. Existence and local uniqueness

Under the previous assumptions, we prove the following result.

Theorem
( 1 )(Existence
( and
) Local Uniqueness for extended McKendrick Equation with Integral Formulation). Let 𝜙0 be an initial condition
in 𝐿1 ∩ 𝐿∞ R+ , 𝐿1 (D) , 𝐾 a global operator satisfying assumptions (5)–(7), and 𝑏 a natality rate( function
) satisfying (8)–(11). There
exists 𝑇 > 0 and 𝑟 > 0 such that if we define 𝑓0 on 𝑈𝑇 by 𝑓0 (𝑡, 𝑎) = 𝜙0 (𝑎 + 𝑡), then 𝑓0 is in 𝐿∞,1 𝑈𝑇 and there exists a unique 𝑓 in
( )
𝐵 𝑓0 , 𝑟 solution of (4).

Proof (sketch). The proof is inspired by the proof of Cauchy–Lipschitz theorem (see Teschl [27, p.47]), using the Banach fixed-point
theorem (see Agarwal ( et) al. [28]). We denote by 𝐹 (𝑓 ) the mapping between (𝑡, 𝑎) and the right member of (4). We define 𝐹 for 𝑓
on a closed ball 𝐵 𝑓0 , 𝑟 centered on 𝑓0 , with radius 𝑟 judiciously chosen later, using hypotheses on 𝐾 and 𝑏 locally semi-Lipschitz
( ) ( )
character. We firstly prove that 𝐹 is well-defined and that 𝐹 (𝑓 ) is in 𝐵 𝑓0 , 𝑟 for 𝑓 in 𝐵 𝑓0 , 𝑟 .
For the well-definition, we show that
+∞ ( ( ) )
‖𝐹 (𝑓 )(𝑡, 𝑎)‖𝐿1 d𝑎 ⩽ (𝐵 𝑇 + 1) ‖ ‖ ‖ ‖
‖𝜙0 ‖𝐿1 + 𝑐0 ‖𝑓 − 𝑓0 ‖∞,𝐿1 + 𝑀0 𝑇
∫−𝑡
( ) [ ]
where 𝑇 is chosen small enough, 𝑐0 is given ( by)(6) on a neighborhood of 0, 𝜙0 and 𝑀0 is given by (7) on 0, 𝑟0 × {𝜙0 }, where 𝑟0
is the radius of an open ball centered on 0, 𝜙0 included in a such neighborhood.
( )
To show the stability of 𝐵 𝑓0 , 𝑟 by 𝐹 , we firstly show that
( ) ( )
𝑇 𝐵‖ ‖ ‖ ‖
‖𝜙0 ‖𝐿1 + ‖𝜙0 ‖𝐿∞ + 𝑇 𝑐0 𝑟 + 𝑀0 + 𝜀0
‖𝐹 (𝑓 ) − 𝑓0 ‖ 1 ⩽
‖ ‖∞,𝐿 1 − 𝐵𝑇
𝑟0
where 𝜀0 is chosen as . We suppose that 𝑟 and 𝑇 are chosen such that
( 4 )
𝐵‖ ‖ ‖ ‖
‖𝜙0 ‖𝐿1 + ‖𝜙0 ‖𝐿∞ + 𝑀0 𝑇 + 𝜀0
( ) ⩽ 𝑟,
1 − 𝐵 + 𝑐0 𝑇
which gives
‖𝐹 (𝑓 ) − 𝑓0 ‖ 1 ⩽ 𝑟.
‖ ‖∞,𝐿
( )
Thus, 𝐵 𝑓0 , 𝑟 is stable by 𝐹 .
( ) ( )
We then show that 𝐹 is contractant on 𝐵 𝑓0 , 𝑟 . More precisely, we show that for 𝑓 and 𝑔 in 𝐵 𝑓0 , 𝑟 ,
( )
( ) ( ( )) 𝑇 2
‖𝐹 (𝑓 ) − 𝐹 (𝑔)‖∞,𝐿1 ⩽ 𝑐0 + 𝑐1 ‖ ‖
‖𝜙0 ‖𝐿1 𝑇 + 𝐵 𝑐0 + 𝑐1 𝑐0 𝑟 + 𝑀0 ‖𝑓 − 𝑔‖∞,𝐿1 ,
2
( )
where 𝑐1 is given by (10) on a neighborhood of 0, 𝜙0 . We suppose that 𝑟 and 𝑇 are chosen such that
( ) ( ( )) 𝑇 2
𝑐0 + 𝑐1 ‖ ‖
‖𝜙0 ‖𝐿1 𝑇 + 𝐵 𝑐0 + 𝑐1 𝑐0 𝑟 + 𝑀0 < 1.
2
( ) ( )
Thus, 𝐵 𝑓0 , 𝑟 is stable by 𝐹 and 𝐹 is contractant on 𝐵 𝑓0 , 𝑟 . Applying the fixed-point Banach theorem, there exists a unique
( )
fixed-point of 𝐹 in 𝐵 𝑓0 , 𝑟 , which is solution of (4) for 𝑡 in [0, 𝑇 ] and almost everywhere in 𝑎. □

118
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 1. Model generator’s working principle.

2. Model generator

2.1. Motivation

When studying structured population dynamics, our goal was to take into account population structuration to reach a more
precise model. For example, because of cultural reasons, the mortality rate may vary according to the geographical area considered,
other things being equal. When modeling an epidemic, factors such as age, gender, cultural framework or educational level influence
the probability to be infected. It is still possible to model populations and interactions as homogeneous blocks, but the results would
be rough, and the conclusions would not be applicable. This justifies the interest for structured population models.
Building and simulating a model is usually made of several steps. The modeling itself consists in specifying a set of choices
on formalism, hypothesis scopes which states are defined in our model, and how the individual parameters influence transitions
between these states. This abstract model is then transcribed into a formal model, for example a set of UML diagrams, which may
itself give an executable model written in a modeling or programming language.
A modeling language such as Modelica [25] enables users to describe a model as a list of algebraic–differential equations, without
the need to go through translation and simulation steps, where the model description is transformed into numerical results through
the execution of powerful solvers. Thus, users do not need to handle discretization and integration of the ordinary differential
equations, approximate resolutions of algebraic equations, and all the numerical issues that such topics arise. The only input required
by the solver is the model description.
However, writing a model description for large population models, with lots of structuration criteria and values, is less easy
than meets the eye. Avoiding errors (syntax, consistency between equations, model completeness) is tricky. Enhancing a model
description which is already written is painful, time-consuming, and this work is not easily reusable when one wants to multiply
modifications.
Here the need of a model generating tool arises. This tool must be able to take as input a textual and comprehensive description
of a population dynamics model. It will use it to generate model description written in the Modelica language. Then, to change the
structuration hypothesis, for instance, one only has to change the textual input and give it to the generator, which will rewrite the
whole model description corresponding to the new modeling choices.
A model generating tool could be helpful when choosing discretization schemes for the numerical simulation of the model.
Indeed, if we decide to use only an ODE or DAE formulation for the model, we have to apply to (2a)–(2c) a semi-discretization for
the age partial derivative. For example, with an upwind finite difference scheme, it would write:
𝜙 (𝑡, 𝑥) − 𝜙𝑖−1 (𝑡, 𝑥) ( )
d𝑡 𝜙𝑖 (𝑡, 𝑥) + 𝑖 + 𝐾 𝑡, 𝜙(𝑡), 𝑎𝑖 , 𝑥 = 0 (12)
𝛥𝑎
with 𝑖 the number of the age class considered. But the choice of scheme can be seen as arbitrary, and one could want to change the
age discretization scheme without having to rewrite the whole model. A tool that separates the model and the input data from the
numerical scheme could offer a such feature.

2.2. Methods

Our model generator’s working principle is illustrated in Fig. 1. The program is written in the Python language.

An example model. To illustrate how the generator enables to create a Modelica model from simple user input, let us take a dummy
example of a population, age-discretized in five age groups, living in a territory divided in three zones 𝐴, 𝐵, 𝐶. People are born,
get old and die of old age only. The natality rate can depend on age, but not on zone. They may change of zone during their life,
through migration whose intensity depends on age and on attractivity of each zone, which can be arbitrarily fixed for illustration
purpose. For example, the graph represented on Fig. 2 represents the connectivity between zones 𝐴, 𝐵 and 𝐶 in our dummy model.
The numerical parameters for the intensity of these flows are chosen arbitrarily and vary with the age class.

119
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 2. Chosen migration dynamics for the dummy model. The indicated numerical values are nominal values which vary according to age.

Fig. 3. Structuration input files for the dummy model. The value of ‘‘max’’ in ‘‘population’’ means that the age of an individual is truncated at 100 years.

2.2.1. Input data

We write as input a textual description of the model, embedding the model structuration, the description of various operators
used in the model, the correlation between various parameters (for example, there is a correlation of many traits between parents
and their newborn, such as geographic area, income level, etc.), and the numerical values such as natality or mortality rates, initial
distribution of traits in the population.
This information is gathered in a set of JSON files, to facilitate the treatment by the generator.

Population structuration. As illustrated on Fig. 3, we write in a first file (Listing 1) the structuration of population as a set of criteria
and values.4 The age structure will be written in a second file (Listing 2), because age plays a particular role in the equations derived
from the McKendrick extended model. One may also consider many population components, with no exchange between them (such
as in a predator–prey model). Then the structuration of each population component will be written in separate files.

Evolution operators. The input data files also include the description of operators, i.e. what is hidden behind the global operator
evolution 𝐾. For example, in our dummy model, we have mortality as a pointwise operator, migration as an exchange operator,
a combination of a kernel and a pointwise operator, and natality which also has a special role in our metamodel. As shown in
Listing 3, we store the names of such phenomena in front of keys describing the type of operator. These key names (‘‘pointwise’’,
‘‘exchange’’, ‘‘boundary’’) refer to various operator classes implemented in the model generator. The values refer to the specific
operator names given in the model.

{
" pointwise " : " mortality " ,
" exchange " : " migration " ,
" boundary " : " natality "
}
Listing 3: Operator description file for the dummy model

An exchange operator, such as migration, needs to be given a connectivity graph to know which values are connected to
each other. We store this information in a separate JSON file (Listing 4). The convention adopted is that keys represent destination
compartments, and the associated values or lists of values represent origin compartments where the population flow may come
from.

4 In the dummy model, there is only one criterion (zone) with three values (𝐴, 𝐵 and 𝐶).

120
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

{
" trait " : " zone " ,
" graph " : {
" A " : [ " B " ],
" B " : [ " A " , " C " ],
"C": ["A"]
}
}
Listing 4: Connectivity graph file for migrations in the dummy model, immediately derived from the graph on Fig. 2

Behavior laws for non-linear operators. A population dynamics model is usually way more complex than our dummy model. Some
relationships between subpopulations may be non-linear. For example, in an epidemic model, the amount of newly infected people
depends on the product between sensible and infected people. Such relationships are written in a separate file, in a mathematical
textual input using natural symbols, or functions defined in another file. Such files are parsed to capture the non-linear relationships
and then translate them in the executable code.

Numerical parameters. Finally, we need to feed our model generator with the numerical values of structured parameters. For
example, with our age and zone structured model, we need a natality values file with the natality rate of each age group, in each
geographic zone (Listing 5). The amount of values needed may explode when the number of structuration criteria increases, thus
one may use a script to generate such values as a combination of statistical data on a population. The initial condition is also given
in a such format.

[
{ " A " : 0.004, " B " : 0.004, " C " : 0.004},
{ " A " : 0.077, " B " : 0.077, " C " : 0.077},
{ " A " : 0.003, " B " : 0.003, " C " : 0.003},
...
]
Listing 5: Numerical values of natality (excerpt)

2.2.2. From input data to virtual model: the semantic modeler

The first step in our model generator is to take the input data and put it together within a single Python object. This virtual
model (see Fig. 4) will have a treelike structure, composed of one or several components which represent separated populations,
which may interact but without exchanging people.
Each component includes operators as instances of the operator class of the same name. Each operator includes structuration
information, behavior laws or diagrams if necessary, and numerical values. The choices of a strongly typed implementation and
an oriented-object development were guided by the will to provide semantics on input data. For example, to call an appropriate
generation method for each type of operator (pointwise, kernel, natality), in terms of equation generations, we needed to develop
specific methods for each type. The virtual model given to the constructor embeds this meaning to ensure that every piece of the
input data will be treated correctly.

<model . operator . pointwiseOperator . LinearAgePointwiseOperator >

name: mortality
structuration : zone : [ ' A ' , ' B ' , ' C ' ]
ageStructuration : max : 100, steps : 5
data: mortality_data
Listing 6: A mortality operator in the virtual model

The Listing 6 shows how a mortality operator is represented in the virtual model. In the dummy model, the mortality operator is a
linear age-structured pointwise operator. Non-linear operators have also been implemented. The structuration for all criteria except
age is represented in the structuration attribute. For some operators, such as exchange operators, user will need to specify the
criteria susceptible to change between source and destination, and the unchanged ones. The ageStructuration attribute carries
information about the maximal age considered (here, 100 years old) and the number of age groups considered (5 in our example).
The data attribute refers to an instance of the Data class, in which the mortality rates are available. This instance is built from
numerical parameters given as input data.
A virtual model can have as many operators as necessary, each operator describing only one phenomenon undergone by the
population. The tree-like structure enables to always have, for the children nodes of an operator, all the information necessary to
provide the next steps of the simulable model generation.

121
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 4. Structure of a virtual model.

Fig. 5. Variables and main equation for the derived dummy model.

2.2.3. From virtual to derived model: the constructor

The second step is to transform our virtual model into a concrete representation of our variables and equations. This step
prepares the generation of Modelica code, as it reproduces in an object-oriented programming context the structure and the objects
of Modelica: a list of variable declarations, and a list of equations involving the declared variables.
The primitive objects are typed to reproduce Modelica’s typing. The structure of the derived model’s equation is strongly linked
to the extended McKendrick formalism.
In Figs. 5–6, we show parts of the derived model we get from our dummy model. The list of the variables involved in the
model is stored as reference (Listing 7) The population variable is implemented as a multi-dimensional array, in order to handle
multi-structuration. Listing 8 corresponds to the main loop of the age semi-discretization :
𝜙 (𝑡, 𝑧) − 𝜙𝑖−1 (𝑡, 𝑧) ( )
d𝑡 𝜙𝑖 (𝑡, 𝑧) + 𝑖 + 𝐾 𝑡, 𝜙(𝑡), 𝑎𝑖 , 𝑧 = 0 (13)
𝛥𝑎
where 𝑧 is the zone considered, and 𝑖 the number of age group. This age semi-discretization has been performed through an upwind
finite difference method, but any other numeric scheme could be used here. There is a strong isomorphism between (13) and Listing
7, knowing that 𝛥𝑎 is here represented by the parameter dAge, in a context of uniform subdivision of the age interval.
𝐾 will be the sum of 𝐾𝜇 and 𝐾𝑧 operators, respectively standing for mortality and migration flows.
The equation in Listing 9 refers to the description of a mortality operator equation, as a pointwise operator:
( ) ( )
𝐾𝜇 𝑡, 𝜙(𝑡), 𝑎𝑖 , 𝑧 = 𝜇 𝑡, 𝑎𝑖 , 𝑧 𝜙𝑖 (𝑡, 𝑧)

The equation in Listing 10 refers to the description of a natality operator equation, with a sum over all age components represented
by +_𝚊𝚐𝚎. The natality variable on the left-hand side is used for the equation of the youngest age group :

122
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 6. Equations describing operators of mortality, natality and migration for the dummy model (excerpts).

∑
𝑁
( )
𝜙0 (𝑡, 𝑧) = 𝑏 𝑡, 𝑎𝑖 , 𝑧 𝜙𝑖 (𝑡, 𝑧).
𝑖=1

Finally, the equation in Listing 11 refers to the description of a migration operator equation. The oldZone attributes corresponds to
the origin zone, and the zone one refers to the destination. The list of terms considered in such equation is given by the connectivity
graph (Listing 4), to avoid writing unexisting flows between compartments. A such equation would be written as:
( ) ( )
𝐾𝑧 𝑡, 𝜙(𝑡), 𝑎𝑖 , 𝑧 = 𝐴 = − 𝜅𝑧 𝑡, 𝑎𝑖 , 𝐵 → 𝐴 𝜙𝑖 (𝑡, 𝑧 = 𝐵)
( )
+ 𝜅𝑧 𝑡, 𝑎𝑖 , 𝐴 → 𝐵 𝜙𝑖 (𝑡, 𝑧 = 𝐴)
( )
+ 𝜅𝑧 𝑡, 𝑎𝑖 , 𝐴 → 𝐶 𝜙𝑖 (𝑡, 𝑧 = 𝐴)

The sign of terms in the right-hand side corresponds to whether the term is an incoming or outcoming migration flow.

2.2.4. From derived model to executable model: the generator

The last step consists in generating a Modelica code simulable by any Modelica solver such as OpenModelica [29] or Dymola [26].
To do so, we used a templating approach using the Jinja2 Python library [30] to convert the derived model into a Modelica code
file. The power of the templating approach consists in iterating inside the template while rendering it, to write each equation needed,
and keep the derived model synthetic. The Fig. 7 shows an excerpt from a Jinja template used to generate a part of the model, and
the corresponding generated model.
The conversion from PDE to ODE is at the core of our discretization methodology. Firstly, the transport term with the age partial
derivative in (2a) is approximated by a discretization scheme to give (12) for instance. The scheme used in (12), an upwind finite
difference model, could be replaced by any other discretization scheme provided that it has the right properties of consistency and
stability. The error estimates are straightforward to compute. It would also be possible to follow the variable 𝜙 along one of its
characteristics, a common practice when studying transport equations.
Secondly, we give an ODE formulation to numeric solvers through the high-level model description enabled by the Modelica
language. It allows one to avoid delving into the details of numeric time integration, which is not the focus of our work. The
numeric methods used, such as Runge–Kutta 4 and DASSL, allow approximation error estimates. For further reading on the DASSL
solver, see [31].
One could wonder the interest of using the Modelica language instead of a usual programming language, and if this choice has
an impact on the simulation process. The Modelica language [25] is a high-level modeling language enabling to describe models as
a collection of ODEs and algebraic equations, which leads to a DAE model. The object-oriented approach, the comprehensiveness of
equations with no need to write in details the numerical schemes, the richness of libraries and the active nature of the community
make it a key asset to simulate DAE models. In Modelica, there is a cost associated with the translation step, which converts the
model description code into an executable code written in C language. But this is largely balanced by the amount of time saved
when achieving multiple simulations.

123
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 7. From Jinja template to generated Modelica code.

2.3. Summary

We developed and implemented a model generator to handle generation of executable models conforming to our extended
McKendrick formalism. These models are written using the Modelica language, from synthetic inputs that can be manually written
using the JSON language, respecting a set of rules that have been documented. We tackled the issue of writing structured models,
with a number of variables and equations which may explode as the structuration is refined and gets more and more complex. Our
model generator handles ODE and PDE formalisms, and generates as output models written in Modelica, but its structure is robust
enough to consider developments with other formalisms and other modeling or programming languages.
For now, the input data interface is not very user-friendly. JSON files have many advantages – they are light, parsing libraries
already exist in almost every programming language, and the standard is complete enough to store any type of data we would need.
This format of input data is acceptable for research uses, as it enables a synthetic representation of a population dynamic model.
But is not suitable for various modeling purposes, which need a simpler and more ergonomic interface to indicate the structure of
the model and equations. For this reason, we plan to implement a UML/SysML interface between the user and the model generator
input. This interface should be user-friendly, for example with drag-and-dropping boxes to represent the steps of a population process
(for example in education, economics or health fields).

124
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 8. Time evolution (a) and phase portrait (b) of the classic predator–prey model.

3. Application to a structured Lotka–Volterra model

3.1. Introduction and model description

The Lotka–Volterra model has been widely studied by mathematicians [32]. Its goal is to model a predator–prey dynamics
within an ecosystem, where there may be a competition for resources and an asymmetric role played by preys and predators. This
model is well known by the mathematical community. For this reason, we found interesting taking it as a toy model to show how
our modeling framework may apply and what interesting results can emerge from the introduction of structuration thanks to the
extended McKendrick formalism.

3.1.1. The original Lotka–Volterra model

To begin with, let us recall the original Lotka–Volterra equations:
{
d𝑡 𝑋(𝑡) = 𝛼 𝑋(𝑡) − 𝛽 𝑌 (𝑡)𝑋(𝑡)
(14)
d𝑡 𝑌 (𝑡) = 𝛿 𝑋(𝑡)𝑌 (𝑡) − 𝛾 𝑌 (𝑡)

where 𝑋 (resp. 𝑌 ) is the number of preys (resp. predators), 𝛼 is the prey natality rate, only due to intrinsic reproduction, 𝛽 is
the prey mortality rate, due to predation, 𝛿 is the predator natality rate, due to predation, and 𝛾 is the predator mortality rate, only
from natural origin.
For the sake of completeness, we remind the classic shapes of time-evolution and phase portrait for the usual predator–prey
model (Fig. 8).

3.1.2. Structuration of prey and predator populations

To illustrate our modeling framework and the model generator’s abilities, we decide to give to our populations of preys and
predators a structuration according to many criteria.

Species structuration. We will separate preys and predators in components (see Fig. 4), as they will not necessarily have the same
structuration and no population exchange will happen between these populations. The prey (resp. predator) component will be
denoted by 𝜙𝑋 (resp. 𝜙𝑌 ).

Age structuration. Both populations will be age structured, to take into account the aging process. The time derivatives in (14) will
be replaced by transport operator 𝜕𝑡 + 𝜕𝑎 .

Territory structuration and migrations. We suppose that the land where our preys and predators live is divided into four zones,
numbered from 1 to 4. We do arbitrary hypotheses on connectivity between zones and numerical parameters representing the
proportion of each zone population which will move to another zone – see Fig. 9. The values indicate a nominal exchange rate,
which is then linearly modulated according to the age and eventually the health status of individuals.

125
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 9. Chosen migration dynamics for the structured Lotka–Volterra model.

Fig. 10. SIRD epidemic model.

Epidemic structuration. For illustration purpose, we decide to assume that the prey population is affected by an epidemic, according
to a SIRD model [11]. Susceptible (𝑆) individuals, once infected, move to the compartment of infected (𝐼) individuals, which may
recover (𝑅) from their disease or die (𝐷). This simple model does not take into account reinfection (no transition from 𝑅 to 𝑆 or
𝐼) nor other mortality factors than the infection, such as natural mortality due to aging. The equations of this model are given by
(15). The compartments are illustrated on Fig. 10.
⎧
⎪ d𝑡 𝑆 = −𝑘𝐼 𝑆
⎪ d𝑡 𝐼 = 𝑘𝐼 𝑆 − (𝜆 + 𝜈)𝐼
⎨ d 𝑅 = 𝜆𝐼 (15)
⎪ 𝑡
⎪ d𝑡 𝐷 = 𝜈 𝐼
⎩

𝑘 represents the infectivity rate, 𝜆 the recovering rate, and 𝜈 the lethality rate.
The goal is to include each of these structuration criteria and the behavior laws governing these phenomena, as they are, for
now, isolated from each other, within a unique population dynamics model, to describe a structured predator–prey model. The
modeling framework developed around the extended McKendrick formalism will tackle this issue.

3.2. Rewriting the problem to conform to the extended McKendrick framework

3.2.1. Extended-McKendrick equations

As explained above, the population variable 𝜙 will be split into two separated components 𝜙𝑋 and 𝜙𝑌 . What follows will still
( )
conform to the McKendrick formulation (3a)–(3b), as we may consider 𝜙 as a couple 𝜙𝑋 , 𝜙𝑌 and study both equations separately.
The hyperparameters set will be, for preys, D𝑋 = [[1, 4]] × {S, I, R, D}, and for predators, D𝑌 = [[1, 4]]. As D𝑋 and D𝑌 are finite sets, the
∑4 ∑
integrals d𝑥 will be replaced by discrete summations on elements of D𝑋 or D𝑌 , for example for preys. We assumed
∫D
𝑧=1 𝑠∈{S,I,R,D}

126
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

that an individual whose parents are located in zone 𝑧 at birth is born in zone 𝑧. We get the following equations:
⎧ 𝜕 𝜙 (𝑡, 𝑎, 𝑧, 𝑠) + 𝜕 𝜙 (𝑡, 𝑎, 𝑧, 𝑠) + 𝐾 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, 𝑠) = 0 ∀ (𝑡, 𝑎, (𝑧, 𝑠)) ∈ 𝐼 × R+ × D𝑋 , (a)
⎪ 𝑡 𝑋 𝑎 𝑋 𝑋
⎪
+∞ ∑
⎨ 𝜙𝑋 (𝑡, 0, 𝑧, 𝑠) = ∫ 𝑏𝑋 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, 𝑠, 𝑠) ̃ d𝑎
̃ 𝜙𝑋 (𝑡, 𝑎, 𝑧, 𝑠) ∀ (𝑡, (𝑧, 𝑠)) ∈ 𝐼 × D𝑋 , (b) (16)
⎪ 0 𝑠∈{S,I,R,D}
̃
⎪ 𝜙𝑋 (0, 𝑎, 𝑧, 𝑠) = 𝜙0,𝑋 (𝑎, 𝑧, 𝑠) ∀ (𝑎, (𝑧, 𝑠)) ∈ R+ × D𝑋 . (c)
⎩
⎧ 𝜕𝑡 𝜙𝑌 (𝑡, 𝑎, 𝑧) + 𝜕𝑎 𝜙𝑌 (𝑡, 𝑎, 𝑧) + 𝐾𝑌 (𝑡, 𝜙(𝑡), 𝑎, 𝑧) = 0 ∀ (𝑡, 𝑧) ∈ 𝐼 × R+ × D𝑌 , (a)
⎪ +∞
⎨ 𝜙𝑌 (𝑡, 0, 𝑧) = ∫ 𝑏𝑌 (𝑡, 𝜙(𝑡), 𝑎, 𝑧) 𝜙𝑌 (𝑡, 𝑎, 𝑧) d𝑎 ∀ (𝑡, 𝑧) ∈ 𝐼 × D𝑌 , (b) (17)
⎪ 0
⎩ 𝜙𝑌 (0, 𝑎, 𝑧) = 𝜙0,𝑌 (𝑎, 𝑧) ∀ (𝑎, 𝑧) ∈ R+ × D𝑌 . (c)

3.2.2. Global operators

For each species, we decompose the global operators 𝐾𝑋 and 𝐾𝑌 as a sum of operators corresponding to each phenomenon,
except for natality which is handled in the boundary condition. For preys, we take into account a 𝐾𝑋 ,𝑚 operator for migrations,
a 𝐾𝑋 ,𝜇 operator for mortality, and a 𝐾𝑋 ,𝑠 operator for epidemic. For predators, we only take into account a 𝐾𝑌 ,𝑚 operator for
migrations and a 𝐾𝑌 ,𝜇 operator for mortality.
{
𝐾𝑋 = 𝐾𝑋 ,𝑚 + 𝐾𝑋 ,𝜇 + 𝐾𝑋 ,𝑠 ,
𝐾𝑌 = 𝐾𝑌 ,𝑚 + 𝐾𝑌 ,𝜇 .

Migration. Let us denote by 𝑚𝑖𝑗 the migration flow intensity from zone 𝑖 to zone 𝑗, as introduced in Fig. 9. The migration operator
can be written as
∑
4
( )
𝐾𝑋 ,𝑚 (𝑡, 𝜙(𝑡), 𝑎, 𝑖, 𝑠) = 𝑚𝑖𝑗 𝜙𝑋 (𝑡, 𝑎, 𝑖, 𝑠) − 𝑚𝑗 𝑖 𝜙𝑋 (𝑡, 𝑎, 𝑗 , 𝑠)
𝑗=1

for preys and similarly as

∑
4
( )
𝐾𝑌 ,𝑚 (𝑡, 𝜙(𝑡), 𝑎, 𝑖) = 𝑚𝑖𝑗 𝜙𝑌 (𝑡, 𝑎, 𝑖) − 𝑚𝑗 𝑖 𝜙𝑌 (𝑡, 𝑎, 𝑗)
𝑗=1

for predators. We observe that each operator does not depend on the other species, and that age and health are supposed not to
change when an individual migrates, as this change can be considered as instant. Such exchange operators can be rewritten as the
sum of a pointwise and of a kernel operator. For example, for predators, we may write
( 4 )
∑ ∑
4
𝐾𝑌 ,𝑚 (𝑡, 𝜙(𝑡), 𝑎, 𝑖) = 𝑚𝑖𝑗 𝜙𝑌 (𝑡, 𝑎, 𝑖) − 𝑚𝑗 𝑖 𝜙𝑌 (𝑡, 𝑎, 𝑗)
𝑗=1 𝑗=1

∑
4
where 𝑚𝑖𝑗 is the rate of the pointwise operator representing the emigration flow from zone 𝑖, and 𝑚𝑗 𝑖 is the value of the (discrete)
𝑗=1
kernel of the kernel operator representing the immigration flow from zone 𝑗 to zone 𝑖.

Mortality. The mortality operators 𝐾𝑋 ,𝜇 and 𝐾𝑌 ,𝜇 will typically be defined as pointwise operators, according to the following
expressions issued from mortality terms in (14):
( +∞ )
𝐾𝑋 ,𝜇 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, 𝑠) = 𝛽 (𝑎, ̃ 𝑧) d𝑎̃ 𝜙𝑋 (𝑡, 𝑎, 𝑧, 𝑠),
̃ 𝑧) 𝜙𝑌 (𝑡, 𝑎,
∫0

𝐾𝑌 ,𝜇 (𝑡, 𝜙(𝑡), 𝑎, 𝑧) = 𝛾 (𝑎, 𝑧) 𝜙𝑌 (𝑎, 𝑧).

The non-dependence of 𝛽 and 𝛾 on 𝑡 is purely a modeling choice.

Epidemic status evolution. The SIRD equations are rewritten in the McKendrick formalism, using pointwise and kernel operators,
resulting in a sum on health statuses.
The sign is reversed in (18a)–(18d) because the global operator 𝐾𝑋 is on the left side of the extended McKendrick equation (16a).
⎧ ( +∞ )
⎪ 𝐾𝑋 ,𝑠 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, S) = 𝑘 (𝑎, 𝑎) ̃ 𝑧, I) d𝑎̃ 𝜙𝑋 (𝑡, 𝑎, 𝑧, S)
̃ 𝜙𝑋 (𝑡, 𝑎, (a)
⎪ ∫0
( )
+∞
⎪
⎨ 𝐾𝑋 ,𝑠 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, I) = − ∫ 𝑘 (𝑎, 𝑎) ̃ 𝑧, I) d𝑎̃ 𝜙𝑋 (𝑡, 𝑎, 𝑧, S) + (𝜆(𝑎) + 𝜈(𝑎)) 𝜙𝑋 (𝑡, 𝑎, 𝑧, I) (b)
̃ 𝜙𝑋 (𝑡, 𝑎, (18)
⎪ 0
⎪ 𝐾𝑋 ,𝑠 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, R) = −𝜆(𝑎)𝜙𝑋 (𝑡, 𝑎, 𝑧, I) (c)
⎪ 𝐾𝑋 ,𝑠 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, D) = −𝜈(𝑎)𝜙𝑋 (𝑡, 𝑎, 𝑧, I) (d)
⎩

In (18a) the flow of susceptible individuals is represented as a pointwise operator, with a rate expressed as a continuous sum
over ages of infected individuals a susceptible individual could meet. 𝑘 (𝑎, 𝑎)
̃ is the infectiosity between infected individuals at age

127
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

𝑎̃ and susceptible individuals at age 𝑎. In (18b) the flow of infected individuals is the sum of two terms. The first one is a kernel
operator representing the incoming infection flow due to contacts between susceptible and infected individuals, same as in (18a)
but with opposite sign. The second one is a pointwise operator representing the outgoing flows of infected individuals who recover,
with recovery rate 𝜆(𝑎), or die, with lethality 𝜈(𝑎). The Eqs. (18c) and (18d) are kernel operators for incoming flows of recovered
or died infected individuals.
Behind such equations, many modeling hypotheses have been done. For example, infectiosity, recovery rate and lethality do not
depend on time 𝑡, nor on geographic zone 𝑧. We observe that in kernel operators, kernel can be very sparse, because the incoming
flow may come from a very restricted sample of compartments. In (18c)–(18d), the incoming flow comes only from the infected
compartment of the same zone, resulting in a very sparse kernel.

3.2.3. Natality
Let us develop the natality rates behind 𝑏𝑋 and 𝑏𝑌 in (16b) and (17b).

Prey natality rate. We suppose that the prey natality rate is constant over time, dependent on the age 𝑎 of the parent, and independent
of the global population state 𝜙(𝑡) and of the zone of birth 𝑧. It may depend on the parent’s epidemic status 𝑠,̃ as infected individuals
may be less fertile, and individuals who died from the epidemic cannot reproduce. We suppose that individuals are born in the S
compartment only. There exists a mapping 𝛼 such as
̃ = 𝛿𝑠S 𝛼 (𝑎, 𝑠)
𝑏𝑋 (𝑡, 𝜙(𝑡), 𝑎, 𝑧, 𝑠, 𝑠) ̃

Predator natality rate. The formulation of the predator natality rate requires to take into account the non-linearity of the dependence
on prey populations. We will then have
+∞ ∑
𝑏𝑌 (𝑡, 𝜙(𝑡), 𝑎, 𝑧) = 𝛿 (𝑎, 𝑎, ̃ 𝑧, 𝑠) d𝑎.
̃ 𝑧, 𝑠) 𝜙𝑋 (𝑡, 𝑎, ̃ (19)
∫0
𝑠∈{S,I,R,D}

(19) expresses the contribution of preys of any age 𝑎̃ and epidemic status 𝑠, in the fixed and common zone 𝑧, to the natality rate
of predator with age 𝑎. This is a typical case where the global population state 𝜙(𝑡) at instant 𝑡 is necessary to compute the natality
rate, through the prey population 𝜙𝑋 (𝑡) specified in an age 𝑎̃ and a hyperparameter (𝑧, 𝑠).

3.3. Numerical results

The following results aim to show the diversity of behaviors one can encounter in a structured population dynamics model,
and the difficulty in forecasting a dynamic system’s asymptotic behavior. We base this statement on three numerical studies on
structured Lotka–Volterra models. The first model consists in a set of 400 variables and equations, representing preys and predators
with 20 age classes on each species, 4 zones, and 4 health statuses for preys. In the second model we removed the epidemic aspect,
took 30 age classes and kept geographical structuration parameters same. The model results in 240 variables and equations. In the
third model, we did not consider any structuration beyond the number of species, so that we only had 4 homogeneous species and
then, 4 variables and equations.
These three examples illustrate how one may numerically grasp the notion of bifurcation in structured population models. Not
only epidemic but also predation parameters, and even age discretization can generate bifurcations. Such bifurcations can be easily
located and visualized using a model generator.
The models generated can be found as supplementary material (Online Resources 1, 2, 3).

3.3.1. When epidemic switches asymptotic behaviors

In the usual predator–prey model studied by Lotka and Volterra, solutions are known to be periodic [33]. Adding spatial
structuration and migration between zones can change the type of asymptotic behavior of the aggregated variables (the total amount
of preys, and the total amount of predators) and present convergence to a fixed point. An epidemic could act as a perturbation of
the system and preserve the fixed point convergence behavior. Surprisingly, it has been observed through numerical investigations
that for a precise parameter set of natality and mortality for preys and predators, the presence of an epidemic could asymptotically
make the system periodic again (see Fig. 11).
For others demographic parameters, we observed that the populations, in a no-epidemic scenario, adopted a periodic trajectory,
whereas the presence of an epidemic changed the behavior to convergence to a fixed point (see Fig. 12).
Two statements may be inferred from these examples, which are not isolated. Firstly, it is very difficult to forecast the future
of a structured population. Many counter-intuitive phenomena can occur, and some parameters can have unexpected influence
on the shapes of solutions. In our examples, we showed how demographic parameters affect both the asymptotic behavior of the
population without the epidemic and the influence of the epidemic on this behavior. Secondly, the power of a model generator
approach is enlightened, with the ability to build structured models over and over again until we find the appropriate structure and
the appropriate parameters to study precise behaviors.

128
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 11. Predator–prey phase portrait. The presence of an epidemic has changed the asymptotic behavior from convergence to periodicity.

Fig. 12. Predator–prey phase portrait. The presence of an epidemic has changed the asymptotic behavior from periodicity to convergence.

3.3.2. Simultaneous existence of stable cycles

Stability of fixed point and limit cycles is major in dynamic systems. It is well known that a dynamical system may have stable
and unstable equilibrium points. The stability can be determined by studying the sign of real parts of the jacobian’s eigenvalues at
the equilibrium point. Equilibrium points partition the space of initial conditions into basins of attraction. Two stable equilibrium
points may exist simultaneously, and the limit of the system will be determined according to which basin of attraction the initial
condition belongs to.
We observed in our structured predator–prey models an interesting, similar phenomenon on stable cycles. The epidemic feature
has been removed from the model for the sake of simplicity. The demographic parameters remained the same, and by changing the

129
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 13. Phase portraits of predator–prey models with different initial conditions (a) and (b). The initial condition has been marked with a red star. The
demographic parameters are unchanged between the two simulations. (For interpretation of the references to color in this figure legend, the reader is referred
to the web version of this article.)

Fig. 14. Time evolution by zone for preys and predators in a prey–predator model.

initial condition, we observed that the system winds around a different cycle (Fig. 13).

3.3.3. Time evolution by zone

Thanks to the model generator, one may plot a detailed evolution of populations over time, selecting one or many criteria. For
example, in a prey–predator model without epidemic dynamic, for some choice of parameters and initial condition, we can visualize
the evolution by species and by zone (see Fig. 14). Both predation dynamic and territorial exchange influence the evolution of the
populations. Such a structured model can provide insights about the detailed evolution and repartition of subpopulations, something
that a homogeneous model could not achieve.

3.3.4. Chaotic behavior

Finally, we highlight a chaotic behavior the predator–prey model may present if one gives it a few complexification. For example,
we consider an ecosystem with no structuration, neither on age nor on geographic zones, but with two prey species (rabbits 𝑥1 and

130
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 15. [llustration of chaos in a two preys–two predators system....]Illustration of chaos in a two preys–two predators system. The bifurcation
diagram (Fig. 15(a)) has been plotted with a 𝜆 parameter varying from 0 to 1 to explore a chosen line segment in the R12 parameter space,
( )
between the following parameter values for 𝛼1 , 𝛼2 , 𝛽11 , 𝛽12 , 𝛽21 , 𝛽22 , 𝛾1 , 𝛾2 , 𝛿11 , 𝛿12 , 𝛿21 , 𝛿22 : (0.5, 0.533, 0.666, 0.03, 0.033, 0.066, 0.046, 0.023, 0.02, 0.046, 5.1, 5.1) and
(0.5, 0.483, 0.042, 0.045, 0.458, 0.042, 0.029, 0.041, 0.035, 0.029, 5.1, 5.1). The Poincaré section (Fig. 15(b)) has been chosen for 𝜆 ≃ 0.208.

squirrels 𝑥2 ) and two predator species (foxes 𝑦1 and wolves 𝑦2 ) living together. The equations of a such system are the following :
⎧
⎪ d𝑡 𝑥1 (𝑡) = 𝛼1 𝑥1 (𝑡) − 𝛽11 𝑦1 (𝑡)𝑥1 (𝑡) − 𝛽12 𝑦2 (𝑡)𝑥1 (𝑡)
⎪ d𝑡 𝑥2 (𝑡) = 𝛼2 𝑥2 (𝑡) − 𝛽21 𝑦1 (𝑡)𝑥2 (𝑡) − 𝛽22 𝑦2 (𝑡)𝑥2 (𝑡)
⎨ (20)
⎪ d𝑡 𝑦1 (𝑡) = 𝛿11 𝑥1 (𝑡)𝑦1 (𝑡) + 𝛿12 𝑥2 (𝑡)𝑦1 (𝑡) − 𝛾1 𝑦1 (𝑡)
⎪ d𝑡 𝑦2 (𝑡) = 𝛿21 𝑥1 (𝑡)𝑦2 (𝑡) + 𝛿22 𝑥2 (𝑡)𝑦2 (𝑡) − 𝛾2 𝑦2 (𝑡)
⎩
where 𝛼𝑖 , 𝛽𝑖𝑗 , 𝛿𝑖𝑗 and 𝛾𝑖 , for 𝑖, 𝑗 = 1, 2, are natality and mortality parameters, defined similarly as in the classic Lotka–Volterra model.
The appearance of chaos in such a simple model may be highlighted through the plot of a bifurcation diagram (Fig. 15(a)). We
let vary parameters in given ranges, and for each value, plotted intersection points of an orbit of (20) with a chosen hyperplane
of R4 . One may distinguish parameter ranges of chaotic behavior, where intersection points (here their projection on the 𝑦1 axis)
spread on a large range of values, from periodic behavior, where intersection points form regular lines. For 𝜆 ≈ 0.8, we observe a
phenomenon that needs explanations. Two periodic attractors coexist, and the initial condition switches from an attraction basin to
the other. On Fig. 15(b), we plotted an example of a Poincaré section. The intersection points with a given hyperplane, for a chosen
set of parameters, were plotted and colored by order of intersection. The various ramifications observed in this figure indicate a
fractal nature that we can relate to usual chaotic systems (see [34]).

Discussion

The developed model generator is able to generate models conforming to a modeling framework introduced to systematically
embed any population dynamics related problem with a stock-flow approach. The modularity of its implementation makes seamless
changing the numerical scheme used to discretize our equations. As we mentioned in Section 2.3, an objective would be to implement
a user–friendlier interface to input model descriptions into the generator, for example using the SysML modeling language. Let us
give an overview of what such an interface could achieve. The SysML language could be used to replace high-level equations by
symbolic links between diagrams. Building a library of architectural patterns or a SysML profile to implement ODEs, pointwise,
kernel and exchange operators would be interesting, and should embed the behavior laws of each component. Thus, the diagrams
will convey the same information as the JSON files, with an interface where a box can be replaced by a more complex one to refine
the model without affecting the whole structure. The generator should be enhanced to convert a diagram into the corresponding
list of structured equations, according to a structuration setting that would also be formulated as a diagram in the SysML language.
It should be noted that a SysML extension named SysPhS [35] exists and enables conversion of SysML models to Modelica code. It
has been used in a first approach but the limitations encountered following our goal of generating structured model in Modelica led
to the development of our own modeler. The Fig. 16 shows an example of structure for a predator–prey model designed in SysML
to be then converted to Modelica through the model generator.
Using a model generator like the one presented could be enlightening for convergence studies on age-structured models. One
could generate models with a chosen number of age classes and observe the influence of this choice over the numerical results.
It could also be used to put to the test numerical studies on age-structured models, to ensure the number of age classes do not
influence the results.
Another limit is of technical nature. Simulating large models written in Modelica is tough. When tackling the issue of modeling
a complete population, with many levels of structuration, one can easily reach the amount of 100 000 variables and equations to
describe a such model. Generating it does not present any problem, but simulating it through Modelica solvers, as we tried with the
one embedded in the Dymola software, may present crippling computation durations. This is a major obstacle to the calibration of

131
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Fig. 16. Example of what a SysML interface for modeling system dynamics model could look like.

these models, as calibration needs hundreds of simulation runs to converge to an optimal distribution of parameters.
Concerning the numerical results presented in Section 3.3, the structured Lotka–Volterra model, which has been built, simulated
and studied, is a toy model. It aims to demonstrate our model generator’s abilities and how complex behaviors may emerge from a
structured model. Our goal is not to forecast how prey and predator populations would evolve in an ecosystem if the structure was
precisely taken into account. Thus, the model was not calibrated using field data, and many numerical values of parameters in it
were chosen arbitrarily, for demonstration purposes. This toy model has been considered to illustrate our modeler abilities, and we
may easily make it more complete and more complex by refining the modeling, taking into account new phenomena, etc. Again,
this refinement would be helped by the large scope of our modeling framework.
The epidemiologic model could be refined, by differentiating predation coefficients for infected preys, or by integrating a
transmission of the disease to predators. The goal with Lotka–Volterra and SIRD models was to show how usual population dynamics
models can be naturally embedded in our framework.

Future work

We plan several in-depth studies to carry on these works. As mentioned before, a SysML interface for the model generator input
is planned for implementation; this will increase the efficiency in building new models and imagining new scenarios. We also plan
to enhance the model generator to build and simulate models in various modeling paradigms, such as multi-agent system, stochastic
simulation, and even build hybrid models mixing such paradigms on different scales.
Deepen the structured Lotka–Volterra model analysis could be of major interest, and deeply relatable with recent works [36–38],
for example by studying and understanding the influence of an age structure on the modeling of predation between individuals. The
stability study in a structured Lotka–Volterra model could be part of this analysis, to better understand the coexistence of stable
cycles and the bifurcation phenomena observed when switching on and off the epidemic operator.

Conclusion

In this work, we developed an agile population modeling framework based on an extension of the well-known McKendrick
equation, which tackles the introduction of structure in population models. This framework has proven to be sufficiently wide to
embed complex structured models and many classic problems in population dynamics that can be easily coupled together (such as
predator–prey and epidemiologic models). Here is the main interest of our approach. It can be relatively easily extended to other
population or biological models, for example working on size structure in cell division models (see Perthame [24], Michel [39]).
We showed that the variety of phenomena usually considered can be reduced to two types of operators, called pointwise and
kernel operators. The theorem of local existence and uniqueness is based on a fixed-point argument, inspired by the proof of Cauchy–
Lipschitz theorem. We implemented a model generator in Python language, in order to generate various models easily and focus
on simulation results instead of writing the model’sexte code manually. We illustrated the power of this approach on a structured
predator–prey model, with age structuration, territorial dynamics and an epidemic through an SIRD model. The introduction of
structure in such model reveals a variety of behaviors, in terms of asymptotic behavior, stability and appearance of chaos.

132
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

Funding

This work was funded by Dassault Systèmes and the Association Nationale de la Recherche et de la Technologie (ANRT).

CRediT authorship contribution statement

Laurent Attias: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Conceptual-
ization. Vincent Siess: Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Investigation,
Formal analysis, Conceptualization. Stéphane Labbé: Writing – review & editing, Validation, Supervision, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared
to influence the work reported in this paper.

Acknowledgments

This work has been carried out under a collaboration between Sorbonne Université and Dassault Systèmes.

Appendix A. Supplementary information

We provide as supplementary materials the models generated to explore numerical behaviors on predator–prey structured
models. These Modelica files contain the full description of these three models, including the numerical parameters used for
demonstration purpose. More precisely, supplementary_material_1.mo refers to the model used in Section 3.3.1, supple-
mentary_material_2.mo refers to the model used in Section 3.3.2, and supplementary_material_3.mo refers to the
model used in Section 3.3.4.

Data availability

As supplementary material.

References

[1] T. Malthus, An Essay on the Principle of Population, J. Johnson, London, 1798.

[2] P. Magal, S. Ruan (Eds.), Structured Population Models in Biology and Epidemiology, in: Lecture Notes in Mathematics, vol. 1936, Springer Berlin
Heidelberg, Berlin, Heidelberg, 2008, [Link]
[3] O. Arino, A survey of structured cell population dynamics, Acta. Biotheor. 43 (1–2) (1995) 3–25, [Link]
[4] D. Mollison, The structure of epidemic models, in: Epidemic Models: Their Structure and Relation To Data, in: Publications of the Newton Institute,
Cambridge University Press, Cambridge, 1995.
[5] Shripad Tuljapurkar, Ulrich K. Steiner, Steven Hecht Orzack, Dynamic heterogeneity in life histories, Ecol. Lett. 12 (1) (2009) 93–106, [Link]
10.1111/j.1461-0248.2008.01262.x.
[6] Ulrich K. Steiner, Shripad Tuljapurkar, Tim Coulson, Carol Horvitz, Trading stages: Life expectancies in structured populations, Exp. Geront. 47 (10) (2012)
773–781, [Link]
[7] S. Kumagai, M.K. Uyenoyama, Genealogical histories in structured populations, Theor. Popul. Biol. 102 (2015) 3–15, [Link]
01.003.
[8] Marcy K. Uyenoyama, Naoki Takebayashi, Seiji Kumagai, Inductive determination of allele frequency spectrum probabilities in structured populations,
Theor. Popul. Biol. 129 (2019) 148–159, [Link]
[9] J.D. Murray, Mathematical Biology: I. An Introduction, third ed., in: Interdisciplinary Applied Mathematics, vol. 17, Springer New York, NY, US, 2002.
[10] A.G. McKendrick, Applications of mathematics to medical problems, P. Edinb. Math. Soc. 44 (1925) 98–130, [Link]
S0013091500034428.
[11] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. Lond. A 115 (772) (1927) 700–721, http:
//[Link]/10.1098/rspa.1927.0118.
[12] W. Feller, On the integral equation of renewal theory, Ann. Math. Stat. 12 (3) (1941) 243–267, [Link]
[13] R. Bellman, K.L. Cooke, Differential-Difference Equations, The RAND Corporation, New York, 1963.
[14] M.E. Gurtin, R.C. Maccamy, Non-linear age-dependent population dynamics, Arch. Ration. Mech. Anal. 54 (3) (1974) 281–300, [Link]
BF00250793.
[15] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, in: Applied mathematics monographs, (no. 7) Giardini editori e stampatori, Pisa,
1995.
[16] O. Diekmann, P. Getto, Boundedness, global existence and continuous dependence for nonlinear dynamical systems describing physiologically structured
populations, J. Differential Equations 215 (2) (2005) 268–319, [Link]
[17] O. Diekmann, et al., The ‘‘cumulative’’ formulation of (physiologically) structured population models, in: P. Clément, G. Lumer (Eds.), Evolution Equations,
Control Theory, and Biomathematics: Proceedings of the Han-sur-Lesse Conference, in: Lecture notes in pure and applied mathematics, (no. 155) Marcel
Dekker, New York, 1994.
[18] O. Diekmann, et al., On the formulation and analysis of general deterministic structured population models I. Linear theory, J. Math. Biol. 36 (4) (1998)
349–388, [Link]
[19] P. Jagers, The Markov structure of population growth, Acta Appl. Math. 14 (1–2) (1989) 103–114, [Link]
[20] P. Jagers, The growth and stabilization of populations, Statist. Sci. 6 (3) (1991) [Link]

133
L. Attias et al. Mathematics and Computers in Simulation 234 (2025) 113–134

[21] O. Diekmann, et al., On the formulation and analysis of general deterministic structured population models II. Nonlinear theory, J. Math. Biol. 43 (2)
(2001) 157–189, [Link]
[22] O. Diekmann, et al., Steady-state analysis of structured population models, Theor. Popul. Biol. 63 (4) (2003) 309–338.
[23] B.L. Keyfitz, N. Keyfitz, The McKendrick partial differential equation and its uses in epidemiology and population study, Math. Comput. Modelling 26 (6)
(1997) [Link]
[24] B. Perthame, Transport Equations in Biology, in: Frontiers in Mathematics, Birkhäuser Basel, Basel, 2007, [Link]
[25] P. Fritzson, V. Engelson, Modelica — A unified object-oriented language for system modeling and simulation, in: E. Jul (Ed.), ECOOP’98 — Object-Oriented
Programming, Springer Berlin Heidelberg, Berlin, Heidelberg, 1998, pp. 67–90, [Link]
[26] D. Brück, et al., Dymola for multi-engineering modeling and simulation, in: M. Otter (Ed.), Proceedings of the 2nd International Modelica Conference,
2002, 55–1–55–8.
[27] G. Teschl, Ordinary Differential Equations and Dynamical Systems, in: Graduate studies in mathematics, (no. 140) American Mathematical Society, RI, US,
2012.
[28] P. Agarwal, M. Jleli, B. Samet, Banach contraction principle and applications, in: Fixed Point Theory in Metric Spaces, Springer Singapore, Singapore,
2018, pp. 1–23.
[29] P. Fritzson, et al., The OpenModelica integrated environment for modeling, simulation, and model-based development, MIC 41 (4) (2020) 241–295,
[Link]
[30] A. Ronacher, Jinja2 Documentation, Technical Report, 2014, URL [Link]
[31] Linda Petzold, A description of DASSL: A differential/algebraic system solver, in: Proceedings of 10th IMACS World Congress, 1982.
[32] P.J. Wangersky, Lotka-Volterra population models, Annu. Rev. Ecol. Syst. 9 (1) (1978) 189–218.
[33] A.J. Lotka, Elements of Physical Biology, in: universalibrary, Williams and Wilkins Company, Baltimore, 1925.
[34] W. Szemplińska-Stupnicka, Chaos: Bifurcations and Fractals Around Us, in: World Scientific series on nonlinear science, vol. 47, World Scientific, New
Jersey, 2003.
[35] Object Management Group, SysML Extension for Physical Interaction and Signal Flow Simulation, Technical Report, 2021.
[36] Y. Lu, K.A. Pawelek, S. Liu, A stage-structured predator-prey model with predation over juvenile prey, Appl. Math. Comput. 297 (2017) 115–130,
[Link]
[37] S. Bentout, S. Djilali, A. Atangana, Bifurcation analysis of an age-structured prey–predator model with infection developed in prey, Math. Methods Appl.
Sci. 45 (3) (2022) 1189–1208, [Link]
[38] X. Zhang, Z. Liu, Hopf bifurcation analysis in a predator-prey model with predator-age structure and predator-prey reaction time delay, Appl. Math. Model.
91 (2021) 530–548, [Link]
[39] P. Michel, Existence of a solution to the cell division eigenproblem, Math. Models Methods Appl. Sci. 16 (2006) 1125–1153, [Link]
S0218202506001480.

134