- 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 denote the concentration of circulating neutrophils in state (with ) at time, where is the number of states identified from the clustering analysis of individual cells. The influx of precursor neutrophils (or precursor state), , is assumed to occur with constant rate , and neutrophils in any state are lost from circulation at rate . A neutrophil transitions from state to state with rate . The changes in neutrophil populations associated with each state may then be described by the following system of ordinary differential equations (ODEs):
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. 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 to be the fraction of circulating neutrophils in state ; that is, . 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
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 be a matrix of these observed fractions, where is the number of subjects we have samples for and is the number of neutrophil states. The row, , then contains the cluster fractions for individual (). For a particular value of , the likelihood of our observed fractions is given by
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 ≤ i ≤ N, 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. 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. 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 (1/ξi,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. 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. 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. 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
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
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
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
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
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
7. Boer RJDe and Yates AJ. Modeling T-cell fate. Annu Rev Immunol. (2023) 41:513–32. doi: 10.1146/annurev-immunol-101721-040924
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
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
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
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
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
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
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
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
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
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 StatesReviewed by:
Isha Monga, Weill Cornell Medicine, United StatesYutaro 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