ABSTRACT (original article): During the twentieth century, inflammatory bowel disease (IBD) was considered a disease of early industrialized regions in North America, Europe and Oceania. At the turn of the twenty-first century, IBD incidence increased in newly industrialized and emerging regions in Africa, Asia and Latin America, while the prevalence in early industrialized regions continued to grow steadily. Changes in the incidence and prevalence denote the evolution of IBD across four epidemiologic stages: stage 1 (emergence), characterized by low incidence and prevalence; stage 2 (acceleration in incidence), marked by rapidly rising incidence and low prevalence; and stage 3 (compounding prevalence), where the incidence decelerates, plateaus or declines while the prevalence steadily increases. A fourth stage (prevalence equilibrium) has been proposed in which the prevalence slope plateaus due to demographic shifts in an ageing IBD population, but it has not yet been evidenced. To date, these stages have remained theoretical, lacking specific numerical indicators to define transition points. Here, using real-world data from 522 population-based studies encompassing 82 global regions and spanning more than a century (1920–2024), we show spatiotemporal transitions across stages 1–3 and model stage 4 progression. Understanding the evolution of IBD across epidemiologic stages enables healthcare systems to better anticipate the future worldwide burden of IBD. CITATION (original article): Hracs, L., Windsor, J.W., Gorospe, J. et al. Global evolution of inflammatory bowel disease across epidemiologic stages. Nature 642, 458–466 (2025). https://doi.org/10.1038/s41586-025-08940-0

Starting Point

This notebook is designed to calculate chronic disease prevalence as a function of time, based on the partial differential equation (PDE) [Brinks et al., ”Age- and time-dependent model of the prevalence of non-communicable diseases and application to dementia in Germany,” Theoretical Population Biology 92, 62-68 (2014)]
∂
∂t
+
∂
∂a
p(t,a)=[1-p(t,a)]i(t,a)-m(t,a)
p(t,a)[R(t,a)-1]
p(t,a)[R(t,a)-1]+1
-p(t,a)r(t,a)+μ(t,a).
(
1
)
Here, a is age, t is time, p(t, a) is prevalence as a function of age and time, i(t, a) incidence rate, m(t, a) population mortality rate, R(t, a) mortality rate ratio (those with disease over those without), r(t, a) recovery rate, and μ(t,a) is migration rate.

Simplified Equation & Solution

We will use a closed-form solution of the prevalence as a function of time and age plus the population age-distribution to create predictions of chronic disease prevalence over time. Here, we can eliminate the recovery term, as one does not typically recover from chronic diseases. To further simplify, we assume that the mortality rate ratio is unity and the migration rate is zero. The simplified PDE for prevalence then takes the form
∂
∂t
+
∂
∂a
p=[1-p(t,a)]i(a,t).
(
2
)
Next, we set
u=1-p
to yield the McKendrick-von Foerster equation
∂
t
u+
∂
a
u=-i(t,a)u,
(
3
)
with 1. initial condition
p(0,a)=
p
0
(a)u(0,a)=1-
p
0
(a)
for some initial age-distribution, and 2. boundary condition
p(t,0)=0u(t,0)=1
, indicating that the prevalence must always go to zero when the age is zero.​Note that for the boundary and initial conditions to be consistent,
p(0,a=0)=
p
0
(0)=0
.​Equation 2 is identical in form to the equation governing age-stratified population changes where the death rate is age- and time-varying and there is a constant (in time and age) birth rate.​Using the method of characteristics, one can solve the equation in two separate regions in age-time space (see page 114, Equation 1 of A.D. Polyanin, V. F. Zaitsev, and A. Moussiaux, “Differential and Integral Equations and Their Applications: Volume 1: Handbook of First-Order Partial Differential Equations,” CRC Press (Taylor & Francis Group), Boca Raton, FL, 2001), allowing both the initial condition and boundary condition to be applied to the solution. ​In the region
0<=t<a
, the initial condition,
u(0,a),
is used to generate the solution
p(t,a)=1-[1-p(0,a-t)]Exp-
t
∫
0
i(x,x+a-t)x,
(
4
)
while in the region
0<=a<t
, the boundary condition,
u(t,0)
, is applied to yield
p(t,a)=1-[1-p(t-a,0)]Exp-
a
∫
0
i(x+t-a,x)x.
(
5
)
The age-independent prevalence at any time can be calculated where there exists age-structured population (i.e., census) data. To do this, we multiply the the prevalence as a function of time and age, p(t, a), by the population age-distribution at time t,
q
t
(a)
, and normalise as per
p(t)=
∞
∫
0
p(t,a)
q
t
(a)a
∞
∫
0
q
t
(a)a
.
(
6
)
The result depends on time only and can be used to extrapolate the prevalence into the future.
​
This model is meant to be conservative when applied to countries with a moderate to high or rapidly growing incidence rate (i.e., the modelled prevalence will likely be greater than what is observed). This is due to two primary factors. First, for countries where where net immigration exceeds net emigration and immigrants typically have a lower chronic disease prevalence than emigrants, the effect will be to slightly over-estimate the prevalence. This is true for most highly industrialised countries. In other regions, where there may be little to no migration and/or the chronic disease burden is similar in the immigrant and emigrant populations, this term is negligible. Second, ignoring the mortality this term will also lead to an over-estimate of prevalence because the mortality rate ratio R(t, a) is greater than unity (higher mortality in the chronic disease population than in the general population).

Required Data

To model the prevalence as a function of time, we need the following data:
1
.
Initial conditions. A CSV containing prevalence at discrete ages must be provided - this serves as the initial condition for the model. The file must have two comma-delimited columns containing the data elements average-age-in-bin and prevalence-at-average-age-in-bin. Since the data are discrete, interpolation is used to calculate the prevalence at any age. To simplify this process, average-age-in-bin should be the mean values between the edges of the age bin. For example, if the age range of the first bin is 0-17, the average-age-in-bin would be 8.5. The code also expects that prevalence-at-average-age-in-bin is a decimal number between 0 and 1.
2
.
Incidence rate. A CSV containing the age-dependent incidence rate for the model. Since there is often much statistical fluctuation in incidence rate data, an average over several years of incidence data could be used. The file must contain two comma-delimited columns, average-age-in-bin and incidence-rate-at-average-age-in-bin. The code expects incidence rates are decimal numbers, not a number per thousand, hundred thousand, etc.
3
.
Population age-structure. One or more CSVs containing the normalised, age-dependent population distribution for the region of interest. For models projecting into the future, these must be estimates of the population age-distribution over the year range of the model. These allow us to integrate out the age dependence in the prevalence to generate a purely time-dependent model. Each CSV must contain two columns: average-age-in-bin and normalised-fraction-at-average-age-in-bin. The normalised-fraction-at-average-age-in-bin value must be a decimal number (the fraction of the population in that age bin normalised to the bin width). Note that the area under the population density curve is normalised to unity in the code.

Loading Data

The rest of the notebook will require the following variables be set.
1) populationYears. Years for which the population age-structure is known.
2) projectionStartYear. The start year for the projection. Currently set to the first year of data in the population distribution data.
3) projectionEndYear. The end year for the projection. Currently set to the last year of the population distribution data set.
4) rootInputDir. The root directory containing the files to load.
5) rootOutputDir. The root directory to which data are to be saved.

Variables

Variables required for loading of the data.
In[]:=
(*Yearsofpopulationdistributiondatatoload.NotethatI'veusedthisvariableinthenamingconventionofthepopulationdistributionfiles.*)​​populationYears=Table[2021+i,{i,0,22}];​​​​(*Startandendyearsforourprojections.*)​​projectionStartYear=populationYears〚1〛;​​projectionEndYear=populationYears〚Length[populationYears]〛;​​​​(*Thisiswhereyouwouldsetinputandoutputdirectories.*)​​​​SetDirectory[NotebookDirectory[]];​​rootInputDir="input/";​​rootOutputDir="output/";

Initial Conditions: Prevalence

Load the age-dependent initial prevalence. Change the prevalencePath variable to select the desired data set.
In[]:=
prevalencePath=rootInputDir<>"prevalence.csv";​​prevalence=Import[prevalencePath,"CSV"];​​prevalence=Table[prevalence〚j〛,{j,2,Length[prevalence]}];​​​​Print["Prevalence: "<>prevalencePath<>" loaded"];​​Print["Number of age data points = "<>ToString[Length[prevalence]]];
Prevalence: input/prevalence.csv loaded
Number of age data points = 11
Alternatively evaluate the following to assign the final form of the data
In[]:=
prevalence=
;

Incidence Rate

Load the age-dependent incidence data. Change the incidencePath variable to select the desired data set.
In[]:=
incidencePath=rootInputDir<>"incidence.csv";​​incidence=Import[incidencePath,"CSV"];​​incidence=Table[incidence〚j〛,{j,2,Length[incidence]}];​​​​Print["Incidence: "<>incidencePath<>" loaded"];​​Print["Number of age data points = "<>ToString[Length[incidence]]];
Incidence: input/incidence.csv loaded
Number of age data points = 11
Alternatively evaluate the following to assign the final form of the data
In[]:=
incidence=
;

Population Age-Distribution

Loading the age-dependent population distribution data. If there is population age-structure data across multiple years, each can be loaded here. In the example, the populationYears variable is used to select the files to be loaded (they are named using the year for which the population distribution what generated).
Set the file path structure by changing the populationPath variable.
In[]:=
populations={};​​For[i=1,i<=Length[populationYears],i++,​​populationPath=rootInputDir<>"population/"<>ToString[populationYears〚i〛]<>".csv";​​population=Import[populationPath,"CSV"];​​population=Table[population〚j〛,{j,2,Length[population]}];​​populations=Append[populations,population];​​​​Print["Population: "<>populationPath<>" loaded"];​​Print["Number of age data points = "<>ToString[Length[population]]];​​];
Alternatively evaluate the following to assign the final form of the data
In[]:=
populations=
;

Interpolations & Plotting

Initial Conditions: Prevalence

Create an interpolating polynomial of the initial prevalence and plot this interpolation. Since we have turned discrete data into continuous data, we need an upper bound for the interpolation, set to the value of the highest age in the given data set.
Set the interpolation order as needed. Be advised that higher order interpolations imply longer model evaluation times and may create unphysical artefacts where the curvature is high.
In[]:=
prevalenceAgeMax=Max[Table[prevalence〚j〛〚1〛,{j,1,Length[prevalence]}]];​​​​interp=Interpolation[prevalence,InterpolationOrder->1];​​prevalenceInterpolation[a_]:=Piecewise[​​{{0,a<0},{interp[a],0<=a<=prevalenceAgeMax},{0,a>prevalenceAgeMax}}​​];​​​​Plot[​​prevalenceInterpolation[a],​​{a,0,prevalenceAgeMax},​​PlotStyle->{Red,Green,Blue},​​Frame->True,ImageSize->450,​​FrameLabel->{Text@Style["Age [Years]"],Text@Style["Prevalence [%]"]},​​LabelStyle->{FontSize->18,FontFamily->"Times",Black},PlotRange->All​​]

Incidence Rate

Create an interpolating polynomial of the incidence and plotting it.

Population Age-Distribution

Create interpolating polynomials of the population age-distribution curves and plot them. Note that the area under the interpolated population distribution is normalised to unity.

Age- & Time-Dependent Prevalence Model

Projected Prevalence

Generate Projections

Calculate the prevalence between the model projectionStartYear and projectionEndYear for one or more incidence growth rates. Set incidenceGrowthRate = {0} to evaluate the model for a constant incidence rate with time only.

Plot Scenarios

Plotting the prevalence as a function of time for the different incidence rate increase or decrease scenarios. Add to the colours variable if there are more than three plots.

Save Prevalence Plots

Save the data used in the prevalence plots in a single file. Change the csvOutputPath variable to adjust the location for the saved plot data.

CITE THIS NOTEBOOK

Global evolution of inflammatory bowel disease across epidemiologic stages​
by Kaplan Global Epidemiology Lab​
Wolfram Community, STAFF PICKS, August 29, 2025
​https://community.wolfram.com/groups/-/m/t/3537476