Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Immunol., 24 October 2025

Sec. Systems Immunology

Volume 16 - 2025 | https://doi.org/10.3389/fimmu.2025.1654015

This article is part of the Research TopicMathematical Modeling in Discovery and Analysis of Immune ResponsesView all 12 articles

A mathematical framework for human neutrophil state transitions inferred from single-cell RNA sequence data

  • 1Systemic Autoimmunity Branch, National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS/National Institutes of Health (NIH)), Bethesda, MD, United States
  • 2Analysis and Intelligence Assessment Directorate, UK Health Security Agency, Salisbury, United Kingdom
  • 3Health Analytics Network, Columbia, MD, United States
  • 4Data Science Unit, The West Bengal National University of Juridical Sciences, Kolkata, India
  • 5Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds, United Kingdom
  • 6Department of Statistics and Applied Probability, University of California Santa Barbara, Santa Barbara, CA, United States
  • 7T-6, Theoretical Biology and Biophysics, Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, United States

Neutrophils, the most abundant immune cells in the human circulation, play a central role in the innate immune system. While neutrophil heterogeneity is a topic of increasing research interest, few efforts have been made to model the dynamics of neutrophil population subsets. We develop a mathematical model to describe the dynamics that characterizes the states and transitions involved in the maturation of human neutrophils. We use single-cell gene expression data to identify five clusters of healthy human neutrophils, and pseudo-time analysis to inform model structure. We find that precursor neutrophils transition into immature neutrophils, which then either transition to an interferon-responsive state or continue to mature through two further states. The key model parameters are the transition rates (the inverse of a transition rate is the mean waiting time in one state before transitioning to another). In this framework, the transition from the precursor to immature state (mean time less than an hour) is more rapid than subsequent transitions (mean times more than 12 hours). Approximately a quarter of neutrophils are estimated to follow the interferon-responsive path; the remainder continue along the standard maturation pathway. We use Bayesian inference to describe the variation, between individuals, in the fraction of cells within each cluster.

Introduction

Neutrophils are the most abundant immune cell subset in the human circulation and play a central role in the rapid innate immune response against microbes and other danger signals (1). Neutrophils are endowed with efficient strategies to neutralize micro-organisms, including phagocytosis, degranulation, and the formation of neutrophil extra-cellular traps (NETs). Conversely, neutrophil dysregulation or enhanced activation can contribute to auto-inflammatory and auto-immune diseases and promote end-organ damage, while decreases in neutrophil function or numbers can lead to enhanced risk for infections (2, 3). The estimated half-life of human mature neutrophils in circulation is considered to be less than a day, which contributes to the challenges of their study and full characterization (4).

Indeed, despite being such a crucial component of the innate immune system, neutrophil physiology and dynamics remain less well characterized than other immune cell types, due in part to technical difficulties in the study of neutrophils, stark differences in neutrophil behavior across species, and the short neutrophil lifetime in circulation. Recent technological advances allow for the investigation of transcriptional states of single neutrophils, which has begun to reveal their heterogeneity and complexity (4, 5). We recently showed that neutrophils can be described along a transcriptional continuum, with discrete clusters based on expression of certain genes (6). However, accurate knowledge of their cell state dynamics remains incomplete, and it is important to characterize the functionality of different neutrophil subsets.

Mathematical modeling has been used to describe the growth, differentiation, and decay of immune cell populations. These models are often based on different hypotheses about the biological mechanisms involved (7). By comparing model predictions with experimental data, researchers cannot only evaluate the plausibility of these hypotheses but also guide the design of future experiments. So far, much of this work has focused on the generation and maintenance of immunological memory, particularly through models of T-cell populations (8). This includes models of CD8+ effector T-cell production, as well as more recent approaches that use single-cell clustering and pseudo-time trajectory analysis to track how tumor cells develop drug resistance (9, 10). In contrast, relatively little effort has been made to mathematically model the transcriptional changes in neutrophils or to quantify their movement through the body beyond what has been observed in vitro (11).

While previous work (6) identified four transcriptionally distinct neutrophil states using single-cell RNA sequencing, it did not provide a quantitative model of their transitions or kinetics. Our study builds on this by introducing a fifth precursor state and developing an integrated mathematical and statistical framework that models the dynamic transitions between neutrophil states. This enables estimation of transition rates, as defined by the mathematical model, and provides a mechanistic interpretation of cell state progression, which was not captured in earlier analyses (6).

Methods

Analysis of scRNA-seq data

Data used for modeling is based on previously published single-cell data from adult human peripheral neutrophils (6). Cells were captured and libraries made using the 10X 3’ 3.1 chemistry. Matrix files were generated using Cellranger 8.0.1 for downstream analysis. The Python package scanpy (version 1.9.4) was used to import matrix files. Cells with fewer than 102 genes or a mitochondrial gene proportion greater than 0.1 were removed. Genes that were expressed in fewer than 5 cells were also omitted from the subsequent analyses.

Healthy volunteer samples were integrated using scvi-tools (version 1.3) and the top 3 × 103 variable genes identified. This creates a low-dimensional, latent representation of the raw gene expression data that attempts to capture the true biological variation while correcting for batch effects. Sample number (1-7) was provided as a categorical covariate to the scVI model, however, donor-specific effects, such as gender, were excluded since their effect cannot be estimated given the small sample size. The latent representation is then used to construct a neighborhood graph that identifies cells that are similar to each other and forms the basis of the clustering. Embedding of this graph in two dimensions is performed using Uniform Manifold Approximation and Projection (UMAP) and allows for better visualization of the clusters. Initial Leiden clustering at a resolution of 0.5 identified eight clusters. Cells in three of these clusters were removed as contaminants due to low expression of CSF3R, FCGR3B and NAMPT, genes that are expected to be highly expressed in human neutrophils. The scVI model was retrained on cells in the remaining clusters and a resolution of 0.42 was used to identify five clusters associated with different neutrophil transcriptional states. A trajectory analysis of these neutrophil states was performed using the Monocle3 Python implementation, py-monocle.

Mathematical model

Let ni(t) denote the concentration of circulating neutrophils in state i (with 1iN) at time, where N is the number of states identified from the clustering analysis of individual cells. The influx of precursor neutrophils (or precursor state), n0, is assumed to occur with constant rate ϕ, and neutrophils in any state are lost from circulation at rate μ. A neutrophil transitions from state i to state j with rate ξi,j. The changes in neutrophil populations associated with each state may then be described by the following system of ordinary differential equations (ODEs):

dn0dt=ϕ(μ+ξ0,1) n0 , dn1dt=ξ0,1n0(μ+ξ1,2+ξ1,3) n1 ,
dn2dt=ξ1,2n1μn2 , dn3dt=ξ1,3n1(μ+ξ3,4) n3 ,
dn4dt=ξ3,4n3μn4.

This system reflects our hypothesis that precursor neutrophils, n0, enter state 1 before differentiating into either state 2 or state 3. The population of neutrophils in state 3 then proceed to the most mature state, represented here by state 4 (see Figure 1).

Figure 1
A UMAP plot on the left displays clusters colored from purple to yellow, labeled 0 to 4, indicating pseudotime progression. A color bar shows pseudotime ranges from 0.0 to 20.0. On the right, a diagram features nodes N0 to N4 connected with arrows labeled with Greek symbols, representing transitions in a process.

Figure 1. Left: UMAP space showing the learned graph and pseudo-time ordering from the trajectory analysis. Right: The mathematical model of the transition of neutrophils between different transcriptional states, derived from the trajectory analysis. Neutrophils enter the circulation from the bone marrow at rate ϕ, with subsequent transitions between states marked by arrows. The state subscripts correspond to the cluster number and the parameter ξi,j represents the rate at which a neutrophil transitions from state Ni to state Nj. Neutrophils in each state are cleared at rate µ, but the arrows depicting these events are omitted from the diagram to improve clarity.

For healthy individuals, we assume that these populations are in equilibrium (stable steady state), and define fi to be the fraction of circulating neutrophils in state i; that is, fi=ni/jnj. The transition rates may then be written in terms of these steady state fractions, all of which are measurable, and the loss rate, μ, as follows

ξ0,1=μ(1f0)f0 , ξ1,2=μf2f1 , ξ1,3=μ(f3+f4)f1 , ξ3,4=μf4f3 .(1)

Note that, for a fixed loss rate, a unique set of cluster fractions yields a unique set of transition rates.

Statistical model and inference

To understand the variability in transition rates across a population of individuals, we begin by using a Dirichlet distribution, Dir(α), to model the population-level variability in steady-state neutrophil fractions. A random variable following a Dirichlet distribution takes values in the form of N-dimensional vectors, with the conditions that each entry of the vector is non-negative and all entries sum to one. Since these conditions apply to the steady state fractions of neutrophils, the Dirichlet distribution provides a natural way to model the between-individual variability in these fractions.

The aim here is, therefore, to estimate the parameters of the Dirichlet distribution, α, using the cluster fractions derived from the analysis of the single-cell sequencing data, our observed fractions. Let F be a M×N matrix of these observed fractions, where M=7 is the number of subjects we have samples for and N=5 is the number of neutrophil states. The mth row, F(m), then contains the cluster fractions for individual m (1mM). For a particular value of α, the likelihood of our observed fractions is given by

L(α;F)=m=1Mg(F(m);α),

where g is the probability density function of the Dirichlet distribution.

The frequentist approach of maximum likelihood estimation looks to find a single α that maximizes the above likelihood, or equivalently, the value of α that makes the observed data most probable. Here, however, we adopt a Bayesian approach that treats our parameters as random variables. In this case, the distribution of α in the presence of our data, or the posterior distribution, is assumed to be proportional to the product of our prior beliefs and the above likelihood (12). Often, as is the case here, precise expressions for the posterior distribution are difficult to obtain, so computational methods such as Markov Chain Monte Carlo (MCMC) are used to instead sample from it. The posterior sampling is performed here using the emcee Python package, version 3.1.4 (13), full details of which are provided in the Supplementary Material. To construct distributions for the transition rates, ξi,j, α are first drawn from the MCMC posterior sample. Then, for each draw, steady-state fractions fi, 1 ≤ iN, are sampled from the corresponding Dir(α) distribution and these sampled fractions are substituted into Equation (1). In this way, we account for both the uncertainty in our estimate of α, and the variability in steady-state fractions across a population of individuals.

Results

Clustering of neutrophils from seven human samples at single-cell level has led to the identification of five clusters corresponding to different transcriptional states. These clusters are depicted in Figure 2, where the UMAP method of dimensionality reduction allows us to visualize the similarity between cells in two dimensions. Full details of the process used to derive these clusters is provided in the Methods section. To interpret the clusters, the average gene expression of cells in each cluster is evaluated for a subset of genes (Figure 3). This subset comprises the top marker genes for each of the clusters found here, as well as marker genes from a previous study of the same dataset (6).

Figure 2
A UMAP plot and a bar chart. The UMAP plot shows colored clusters labeled zero to four, with colors red, yellow, green, blue, and purple. The x-axis is UMAP1 and the y-axis is UMAP2. The bar chart displays the fraction of cells in each cluster by sample HC1 to HC7, with colors corresponding to cluster categories in the legend.

Figure 2. Left: A visualization of the five neutrophil clusters that uses UMAP to produce a two-dimensional representation of the raw gene expression data. Right: For each healthy subject (HC1-HC7), the fraction of neutrophils that belong to each cluster is depicted.

Figure 3
Dot plot showing gene expression across clusters 0 to 4 for genes RPS8 to CSF3R. Circle size indicates fraction of cells in a group, and color intensity shows mean expression level, with a scale from light blue (0.0) to dark blue (1.0).

Figure 3. The average expression of select genes in neutrophils in each of the five clusters (states). These expression levels guide the interpretation of each cluster.

Cluster 0 is characterized by elevated expression of ribosomal protein genes (RPL, RPS), suggesting these are precursor neutrophils that have recently entered circulation from the bone marrow. Cluster 1 shows higher expression of S100 family genes and likely represents immature neutrophils. Cluster 2 defines a distinct state, marked by up-regulation of genes associated with a type I interferon (IFN) response. Clusters 3 and 4 correspond to increasingly mature circulating neutrophils, distinguished by elevated expression of CXCR2 and the long non-coding RNAs, MALAT1 and NEAT1, respectively. This interpretation of neutrophil states aligns with our previous study, which identified four transcriptionally distinct clusters corresponding to clusters 1–4 described here (6), as well a previous study capturing mouse neutrophils from different compartments (14).

To better understand neutrophil dynamics across the identified clusters, we developed a mathematical model based on a system of ODEs, a standard approach for describing how quantities evolve over time. In this model, each variable represents the concentration of neutrophils in one of the five transcriptional states. The concentration within a state can increase due to incoming cells from other states and decrease due to transitions to other states or exit from circulation. As such, the model requires a defined ordering of states.

Pseudo-time analysis was used to order neutrophils along a biological continuum based on gene expression profiles, offering a plausible path through transcriptional states (Figure 1, left). Since precursor neutrophils correspond to the earliest transcriptional state, cells in cluster 0 were used to define the root of the pseudo time analysis. We used this pseudo-time ordering to inform the structure of our model (illustrated in Figure 1, right). In this framework, precursor neutrophils (N0) enter the circulation from the bone marrow and transition into immature neutrophils (N1). From there, cells either pass through an IFN-responsive state (N2) or continue to mature through states N3 and N4. The model’s key parameters are the transition rates, ξi,j, representing the rate at which neutrophils move from state Ni to state Nj. The inverse of these rates (1i,j) corresponds to the average time spent in state Ni before transitioning to Nj. High transition rates imply short dwell times, whereas low rates indicate longer persistence in a state. To estimate these parameters, we assume that neutrophil populations in healthy individuals are at steady-state, meaning the number of cells in each state remains constant over time. Under this assumption, algebraic equations can be derived for the steady-state neutrophil fractions in each cluster, fi, that depend solely on the transition rates and a fixed exit rate from circulation. This exit rate is set based on a presumed 12-hour half-life for circulating neutrophils (15). We then estimate the transition rates of the seven healthy individuals by substituting their observed cell fractions (Figure 2, right) into these steady-state expressions.

Since cluster fractions vary across individuals, the corresponding estimates of transition rates also differ. However, due to the small sample size, these estimates alone do not capture how transition rates vary across the broader population. To address this, we developed a statistical model, described in detail in the Methods section, that provides us with predicted probability distributions for the steady-state cluster fractions (Figure 4). Bayesian inference is used to obtain full posterior distributions for the parameters of this statistical model, while samples from this model act as simulated cluster fractions that are representative of a broader population. These samples are subsequently used to construct probability distributions for the transition rates in the mathematical model. The resulting distributions, provided in Figure 5, account for both the uncertainty in the parametrization of the statistical model and the inherent variability in sampling the cluster fractions. Together, the statistical and mathematical models link individual-level data to population-level estimates of neutrophil state transitions.

Figure 4
Five density plots display data distributions for variables \(f_0\) to \(f_4\). Each plot includes red lines for healthy sample cluster fractions, a blue line for the median, and a gray shading indicating the fifty percent credible interval.

Figure 4. Marginal distributions of the Dirichlet distribution used to model the variability in cluster fractions across a wider population. Each MCMC sample provides a different parametrization of the Dirichlet distribution, and therefore a different density for each cluster. The median across these densities is shown by the solid lines, while the shaded 50% credible intervals indicate regions that half of these densities fall within. On each subplot, the seven red ticks along the x-axes represent the cluster fractions for the seven individuals whose samples are included in the single-cell sequencing analysis.

Figure 5
Four density plots show transition rate densities compared to healthy sample transition rates, indicated by blue curves and red markers. Each plot represents a different logarithmic transition rate variable: logξ₀,₁, logξ₁,₂, logξ₁,₃, and logξ₃,₄, with varying x-axis ranges.

Figure 5. Probability distributions for the transition rates, ξi,j, between state i and state j. Solid lines are model predicted densities, constructed by substituting sampled cluster fractions from the Dirichlet distribution into the mathematical expressions for the transition rates, equation (1). These densities therefore reflect the variability in transition rates across a broader population. Red ticks on the x-axes indicates the values of each transition rate for the seven individuals included in the single-cell sequencing analysis. Transition rates are per hour and logarithmic values are base 10.

Median estimates from this combined framework suggest a rapid transition from precursor to immature neutrophils (ξ0,1 = 4.9 h−1), consistent with the low observed fraction of cells in cluster 0. Transition rates between subsequent states occur on slower but comparable time scales: (ξ1,2 = 0.025 h−1, ξ1,3 = 0.076 h−1, ξ3,4 = 0.055 h−1). From these estimates, the fraction of neutrophils that follow the IFN-responsive path (via cluster 2) can be calculated as: ξ1,2/(ξ1,2 +ξ1,3)=0.24. This suggests that approximately 24% of neutrophils transition into the IFN-responsive state, while the remaining 76% continue along the standard maturation pathway.

Discussion

This study contributes to our understanding of neutrophil heterogeneity by introducing a novel mathematical framework to model transcriptional state dynamics in circulating neutrophils in adult healthy subjects. While mathematical models have been widely applied to study other immune cell populations, particularly T cells, relatively few efforts have addressed the transcriptional progression of neutrophils. Here, we leveraged single cell RNA sequencing data from healthy individuals to identify five distinct transcriptional neutrophil states through clustering analysis. Building on previous work, and based on the expression of key marker genes, these clusters were interpreted to represent a progression from precursor cells to more mature states, with one cluster reflecting a type I IFN-responsive population (6, 14).

By integrating single-cell transcriptomic data with a system of ODEs, we provide a quantitative model for neutrophil state transitions. This model formalizes hypotheses about neutrophil maturation pathways, and enables the estimation of transition rates between transcriptional states. Despite a limited number of healthy subject samples, our statistical approach, based on Bayesian inference, allowed us to incorporate inter-individual variability and generate population-level estimates of state dynamics. This feature is especially valuable for extending insights beyond the specific individuals sampled and provides a framework that can scale with larger datasets of individuals affected by a number of disease states or physiological events (1, 3).

A key insight from the model is the rapid transition from the precursor to immature state (ξ0,1), consistent with the low fraction of neutrophils observed in the precursor cluster. However, the inferred values for ξ0,1 are highly sensitive to small variations in the precursor fraction (f0), which is estimated from a sparse population (≤1% of cells). Since ξ0,1 is approximately equal to µ/f0, small fluctuations in f0 can lead to large changes in ξ0,1. Thus, while the exact value of ξ0,1 may be uncertain, the main conclusion, that this transition occurs rapidly, remains robust. In future studies, it may be reasonable to exclude the precursor cluster and return to the four-cluster model of previous work, particularly if more robust estimates of f0 remain difficult to obtain (6).

We used the value µ = log(2)/12h−1, corresponding to a 12-hour half-life in circulation, that is consistent with established estimates. Since each transition rate in the model scales linearly with µ, uncertainty in µ propagates proportionally to the estimated rates. A + 10% change in µ, for example, would result in a corresponding +10% change in all transition rate estimates. The effect that different neutrophil clearance rates have on the distributions of transition rates is provided in the Supplementary Material. Moreover, the current statistical model assumes that all individuals are sampled from a single population. With a larger and more diverse dataset, it would be possible to incorporate covariates such as age, sex, or clinical status to refine these assumptions and better capture population-level variability.

Our current model assumes that neutrophils enter circulation from the bone marrow at constant rate ϕ. As a result, the ODE model reaches equilibrium and the steady state cluster fractions are constant. A more realistic assumption may be to use a time-dependent, periodic function, ϕ(t), that captures the circadian rhythm of circulating neutrophils (16). In this case, the system of ODEs never reaches equilibrium, with the cluster fractions instead oscillating around the constant steady state values predicted by our current model (see Figure 6). With constant ϕ, the steady state fractions may therefore be interpreted as time-averaged fractions, taken over a 24-hour period. Estimates of the transition rates, ξi,j, are therefore also based on the assumption that cluster fractions from each healthy sample are representative of their daily average. While this acts as a simplification, longitudinal data, describing how cluster fractions for each individual change throughout the day, would be required to parametrize a model with a time-dependent influx that better reflects the true biological system.

Figure 6
Six graphs showing functions over time in hours. Graphs depict \( f_0(t) \) to \( f_4(t) \) and \( \phi(t) \). Blue lines indicate time-dependent \(\phi\), red lines for constant \(\phi\). Functions exhibit exponential decay, oscillation, and steady-state behaviors. Graph \( \phi(t) \) illustrates oscillations, with the red line steady at two billion.

Figure 6. A comparison between the predicted cluster fractions when the influx of neutrophils from the bone marrow, ϕ, is assumed to be constant (red) and varies throughout the day (blue). Here, the time-varying ϕ takes the form of a simple sinusoidal function that describes the circadian rhythm exhibited by circulating neutrophil populations.

Importantly, our model examines and assumes a steady-state neutrophil population (in the different states) for healthy individuals. Extending the framework to capture other non-equilibrium dynamics, such as those occurring in acute infection, chronic inflammation, or drug treatments, represents a valuable future direction. There is, for example, evidence from murine models that the IFN-related neutrophil cluster can expand during bacterial infection, suggesting a functional link in this state (14). Since IFN is involved in both infections and auto-immune conditions, we can speculate that this cell cluster is pre-primed to be reactive and possibly important for disease pathogenesis. Interestingly, it is mostly type I IFN genes that are dynamically regulated in neutrophils (17). In addition, while our model is built on transcriptional data, integrating complementary data types (e.g., proteomics, epigenetics, functional assays, or tissue localization) would provide a more comprehensive view of neutrophil function and state transitions (2, 11). We point out that our cluster interpretation is consistent with previous work (6), yet, the current study goes further by placing these clusters within a dynamic mathematical model. This framework provides not only a snapshot of neutrophil heterogeneity, but also insight into the timing and variability of transitions between states, offering a new dimension to understanding neutrophil maturation.

In summary, this study presents a framework for modeling neutrophil transcriptional dynamics, addressing a significant gap in immunological modeling. Beyond improving our understanding of neutrophil behavior in healthy individuals, this approach lays the groundwork for investigating how neutrophil state dynamics is altered in disease contexts, including autoimmune and auto-inflammatory disorders, infections, metabolic diseases, and cancer. Modeling allows for comparing cell states and transition rates in very different contexts and compartments, to identify commonalities and unique clusters. Hopefully it could guide research and potential treatments targeting specific neutrophil populations in disease. Ultimately, the integration of mathematical modeling with high-dimensional single-cell data can inform both mechanistic understanding and experimental design, advancing the study of neutrophil biology and its role in immune system function.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: GSE188288 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE188288. Python code for the analysis of the single-cell sequencing data and model inference is available at https://github.com/jontycarruthersphe/neutrophil_frontiers_2025.

Author contributions

GW: Conceptualization, Data curation, Investigation, Writing – original draft, Writing – review & editing. JC: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. SR: Data curation, Formal analysis, Methodology, Writing – review & editing. TF: Supervision, Writing – review & editing. GL: Conceptualization, Supervision, Writing – review & editing. SP: Conceptualization, Data curation, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing. CM-P: Conceptualization, Supervision, Writing – original draft, Writing – review & editing. MK: Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. SP was supported by the NIH 2022 National Service DATA Scholarship. The study was also funded in part by the Intramural Research Program at NIAMS/NIH (ZIA AR-041199).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The contributions of the NIH author(s) were made as part of their official duties as NIH federal employees, are in compliance with agency policy requirements, and are considered Works of the United States Government. However, the findings and conclusions presented in this paper are those of the author(s) and do not necessarily reflect the views of the NIH or the U.S. Department of Health and Human Services.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1654015/full#supplementary-material

References

1. Burn GL, Foti A, Marsman G, Patel DF, and Zychlinsky A. The neutrophil. Immunity. (2021) 54:1377–91. doi: 10.1016/j.immuni.2021.06.006

PubMed Abstract | Crossref Full Text | Google Scholar

2. Aroca-Crevillén A, Vicanolo T, Ovadia S, and Hidalgo Andrés. Neutrophils in physiology and pathology. Annu Rev Pathol: Mech Dis. (2024) 19:227–59. doi: 10.1146/annurev-pathmechdis-051222-015009

PubMed Abstract | Crossref Full Text | Google Scholar

3. Wigerblad G and Kaplan MJ. Neutrophil extracellular traps in systemic autoimmune and autoinflammatory diseases. Nat Rev Immunol. (2023) 23:274–88:2023. doi: 10.1038/s41577-022-00787-0

PubMed Abstract | Crossref Full Text | Google Scholar

4. Palomino-Segura M, Sicilia J, Ballesteros Iván, and Hidalgo Andrés. Strategies of neutrophil diversification. Nat Immunol. (2023) 24:575–84. doi: 10.1038/s41590-023-01452-x

PubMed Abstract | Crossref Full Text | Google Scholar

5. Ng LG, Ostuni R, and Hidalgo A. Heterogeneity of neutrophils. Nat Rev Immunol. (2019) 19:255–65. doi: 10.1038/s41577-019-0141-8

PubMed Abstract | Crossref Full Text | Google Scholar

6. Wigerblad G, Cao Q, Brooks S, Naz F, Gadkari M, Jiang K, et al. Single-cell analysis reveals the range of transcriptional states of circulating human neutrophils. J Immunol. (2022) 209:772–82. doi: 10.4049/jimmunol.2200154

PubMed Abstract | Crossref Full Text | Google Scholar

7. Boer RJDe and Yates AJ. Modeling T-cell fate. Annu Rev Immunol. (2023) 41:513–32. doi: 10.1146/annurev-immunol-101721-040924

PubMed Abstract | Crossref Full Text | Google Scholar

8. Antia R, Ganusov VV, and Ahmed R. The role of models in understanding CD8+ T-cell memory. Nat Rev Immunol. (2005) 5:101–11. doi: 10.1038/nri1550

PubMed Abstract | Crossref Full Text | Google Scholar

9. Hamlet Chu H, Chan S-W, Gosling JP, Blanchard N, Tsitsiklis A, Lythe G, et al. Continuous effector CD8+ T cell production in a controlled persistent infection is sustained by a proliferative intermediate population. Immunity. (2016) 45:159–71. doi: 10.1016/j.immuni.2016.06.013

PubMed Abstract | Crossref Full Text | Google Scholar

10. Magi S, Ki S, Ukai M, Domínguez-Hüttinger E, Naito AT, Suzuki Y, et al. A combination approach of pseudotime analysis and mathematical modeling for 281 understanding drug-resistant mechanisms. Sci Rep. (2021) 11:18511. doi: 10.1038/s41598-021-97887-z

PubMed Abstract | Crossref Full Text | Google Scholar

11. Chen W, Boras B, Sung T, Hu W, Spilker ME, and D’Argenio DZ. A whole body circulatory neutrophil model with application to predicting clinical neutropenia from in vitro studies. CPT: Pharmacometrics Syst Pharmacol. (2021) 10:671–83. doi: 10.1002/psp4.12620

PubMed Abstract | Crossref Full Text | Google Scholar

12. Hines KE. A primer on bayesian inference for biophysical systems. Biophys J. (2015) 108:2103–13. doi: 10.1016/j.bpj.2015.03.042

PubMed Abstract | Crossref Full Text | Google Scholar

13. Foreman-Mackey D, Hogg DW, Lang D, and Goodman J. emcee: the mcmc hammer. Publications Astronomical Soc Pacific. (2013) 125:306. doi: 10.1086/670067

Crossref Full Text | Google Scholar

14. Xie X, Shi Q, Wu P, Zhang X, Kambara H, Su J, et al. Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection. Nat Immunol. (2020) 21:1119–33. doi: 10.1038/s41590-020-0736-z

PubMed Abstract | Crossref Full Text | Google Scholar

15. Pérez-Figueroa E, Álvarez-Carrasco P, Ortega E, and Maldonado-Bernal C. Neutrophils: Many ways to die. Front Immunol. (2021) 12:631821. doi: 10.3389/fimmu.2021.631821

PubMed Abstract | Crossref Full Text | Google Scholar

16. Aroca-Crevillén A, Adrover JoséM, and Hidalgo Andrés. Circadian features of neutrophil biology. Front Immunol. (2020) 11:576. doi: 10.3389/fimmu.2020.00576

PubMed Abstract | Crossref Full Text | Google Scholar

17. Grieshaber-Bouyer R, Radtke FA, Cunin P, Stifano G, Levescot Anaïs, Vijaykumar B, et al. The neutrotime transcriptional signature defines a single continuum of neutrophils across biological compartments. Nat Commun. (2021) 12:2856. doi: 10.1038/s41467-021-22973-9

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: neutrophil, cell population, state transition, mathematical model, ODE, single-cell RNA sequence data

Citation: Wigerblad G, Carruthers J, Ray S, Finnie T, Lythe G, Pyne S, Molina-París C and Kaplan MJ (2025) A mathematical framework for human neutrophil state transitions inferred from single-cell RNA sequence data. Front. Immunol. 16:1654015. doi: 10.3389/fimmu.2025.1654015

Received: 25 June 2025; Accepted: 29 September 2025;
Published: 24 October 2025.

Edited by:

Penelope Anne Morel, University of Pittsburgh, United States

Reviewed by:

Isha Monga, Weill Cornell Medicine, United States
Yutaro Kumagai, National Institute of Advanced Industrial Science and Technology (AIST), Japan

Copyright © 2025 Wigerblad, Carruthers, Ray, Finnie, Lythe, Pyne, Molina-París and Kaplan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Saumyadipta Pyne, c3B5bmVAdWNzYi5lZHU=; Carmen Molina-París, bW9saW5hLXBhcmlzQGxhbmwuZ292; Mariana J. Kaplan, bWFyaWFuYS5rYXBsYW5AbmloLmdvdg==

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.