- Agarwala, S., Kim, K.-Y., Phan, S., Ju, S., Ye Eun Kong, Castillon, G. A., Bushong, E. A., Ellisman, M. H., & Tamplin, O. J. (2022). Defining the ultrastructure of the hematopoietic stem cell niche by correlative light and electron microscopy. ELife, 11. https://doi.org/10.7554/elife.64835
- Dingli, D., & Michor, F. (2006). Successful Therapy Must Eradicate Cancer Stem Cells. Stem Cells, 24(12), 2603–2610. https://doi.org/10.1634/stemcells.2006-0136
- Human Recombinant IL-3 | STEMCELL Technologies. (2016). Stemcell.com. https://www.stemcell.com/products/human-recombinant-il-3.html
- Khan, A. O., Rodriguez-Romera, A., Reyat, J. S., Olijnik, A.-A., Colombo, M., Wang, G., Wen, W. X., Sousos, N., Murphy, L. C., Grygielska, B., Perrella, G., Mahony, C. B., Ling, R. E., Elliott, N. E., Karali, C. S., Stone, A. P., Kemble, S., Cutler, E. A., Fielding, A. K., & Croft, A. P. (2023). Human Bone Marrow Organoids for Disease Modeling, Discovery, and Validation of Therapeutic Targets in Hematologic Malignancies. Cancer Discovery, OF1–OF22. https://doi.org/10.1158/2159-8290.cd-22-0199
- Pedersen, R. K., Andersen, M., Stiehl, T., & Ottesen, J. T. (2023). Understanding Hematopoietic Stem Cell Dynamics—Insights from Mathematical Modelling. Current Stem Cell Reports, 9(1), 9–16. https://doi.org/10.1007/s40778-023-00224-5
- Rivera, T., Zhao, Y., Ni, Y., & Wang, J. (2020). Human‐Induced Pluripotent Stem Cell Culture Methods Under cGMP Conditions. Current Protocols in Stem Cell Biology, 54(1). https://doi.org/10.1002/cpsc.117
- Rossi, L., Challen, G. A., Sirin, O., Lin, K. K.-Y., & Goodell, M. A. (2011). Hematopoietic Stem Cell Characterization and Isolation. Methods in Molecular Biology (Clifton, N.J.), 750, 47–59. https://doi.org/10.1007/978-1-61779-145-1_3
- Straube, T., Müller, C., & Perestrelo, A. R. (2023). How to do a Proper Cell Culture Quick Check. Www.leica-Microsystems.com. https://www.leica-microsystems.com/science-lab/life-science/how-to-do-a-proper-cell-culture-quick-check/
- Wojciech Fidyk, Iwona Mitrus, Agnieszka Ciomber, Andrzej Smagur, Agata Chwieduk, Głowala-Kosińska, M., & Giebel, S. (2017). Evaluation of proinflammatory and immunosuppressive cytokines in blood and bone marrow of healthy hematopoietic stem cell donors. Cytokine, 102, 181–186. https://doi.org/10.1016/j.cyto.2017.09.001
- Yadav, P., Vats, R., Bano, A., & Bhardwaj, R. (2020). Hematopoietic Stem Cells Culture, Expansion and Differentiation: An Insight into Variable and Available Media. International Journal of Stem Cells, 13(3), 326–334. https://doi.org/10.15283/ijsc19157
Overview
Developing bone marrow organoids from induced pluripotent stem cells (iPSC) requires a significant amount of financial and temporal capital to culture and maintain. With this in mind, our team worked to create an in silico model that simulated immune response in iPSC-derived bone marrow organoids. This model provides unique insight into the time evolution of bone marrow organoids and aims to further the understanding of the immune system under septic conditions. In addition to furthering the global understanding of sepsis, our model also allows our wet lab team to better identify potential high-yield experiments by predicting which cell populations are significantly altered in the bone marrow. This initial screening helps identify potential experiments at relatively little cost, before having to spend the time and resources necessary to run an experiment on the organoids themselves.
Using a system of differential equations, our model simulates the growth, death, and differentiation of Hematopoietic Stem and Progenitor Cells (HSPCs) into lymphoids, myeloids, and leukocytes. It also models the relative levels of “pro” and “anti” inflammatory biomarkers in response to an infectious injury. As mentioned above, this gives us the ability to better predict the expected trends in cell densities and cytokine concentrations as a result of varying environmental stimuli. Our model differs from other attempts in literature by utilizing an Agent Based Modeling (ABM) approach, a method of programming that focuses on the properties and interactions of specific objects/agents. This allows us to more easily differentiate between different classes of agents (stem cells, differentiated cells, pathogens, cytokines, etc) and assign each class its own unique attributes and characteristics. Using an ABM approach also allows us to output real-time visualizations of the system and its evolution over time, making the model much more accessible to researchers hoping to use the software.
Ideation and Equations
Our model is built on a system of differential equations. These equations are primarily based on existing clinical data, with a few additions and simplifications to tailor the model to fit our needs. Beginning with the basic cell growth equations:
Our model follows 6 different cell populations over time: HSPCs (PH), lymphoids (PL), myeloids (PM), leukocytes (PK), active Leukocytes (PA), and suppressive Leukocytes (PS). HSPCs are the most important population to monitor, as they differentiate into all the other cell types. The growth/death of our HSPC population over time is defined by the HSPC Growth Rate ( GH = 1.005 cells/day), HSPC Death Rate (DH = 0.002 cells/day), and a homeostasis/crowding factor that depends on the current HSPC population in relation to the carrying capacity (Pmax): ɸ = 1/(1+PH⁄Pmax). Using this all together, the rate of change in PH is defined as: dH⁄dt = GHPHɸ - DHPH
The other 5 cell populations were calculated in a different manner. Because they all differentiate from the HSPCs, their populations all depend on the HSPC population and a set differentiation rate. For the lymphoid class of cells, the population PL depends on PH , the HSPC to lymphoid differentiation rate (rHL = 0.005), and the lymphoid death rate DL = 0.15/70. This gives a total rate of change of dL⁄dt = (PH*rHL* ɸ)-(PL DL). The same methodology was used to define the myeloid class of cells. rHM = 0.44/70, DM = 0.15/70 cells/tick → dM⁄dt = (PH*rHM* ɸ)-(PM DM)
The leukocyte class of cells is produced by both lymphoids and myeloids. Lymphoids convert to leukocytes at a rate of rLK=1, meaning every lymphoid converts to a leukocyte. Myeloids convert to leukocytes at a slower rate, approximated to be rMK=0.5. The combined leukocyte growth rate is then defined as the growth rate of PL times the conversion rate rLK plus the growth rate of PM times the conversion rate rMK: → dK⁄dt = dL⁄dt rLK + dM⁄dt rMK
These above equations provide a base for the HSPC culture and differentiation procedures. Next, we defined the immune response equations. The pathogen population (PP) was assumed to grow at a base rate of 1.6 additional pathogens per unit time (GP1). Furthermore, if Pp in a particular cell is greater than 100 (an arbitrary value used for simplicity), an additional pathogen unit is created (GP2). This models the exponential reproduction and spread of pathogens between cells. These relations define the total pathogen growth rate (GP). The pathogen population directly decreases the HSPC population at a rate of [-I * PP ], where I is a scaling constant known as the “infection intensity.” Our model allows this infection intensity variable to be manipulated by users on a scale of 0 to 200. An infection intensity of 0 means there is no pathogenic insult, 100 means a moderate pathogen insult, and 200 is the maximum pathogen insult. Beyond controlling the HSPC death-rate by pathogen, the infection intensity also determines the initial number of pathogens that are introduced to the system. Because of these two factors, the infection intensity scales exponentially.
The growth of pathogen density in a particular cell also impacts the population and differentiation of leukocytes. Leukocytes differentiate into active leukocytes and suppressive leukocytes depending on the pathogen population. Active leukocytes (PA) grow at a rate of dA⁄dt = dK⁄dt * dP⁄dt, and suppressive leukocytes (PS) grow at a rate of dS⁄dt = dK⁄dt ÷ dP⁄dt. The populations of these active and suppressive leukocytes then impact the pathogen population. Active leukocytes reduce the pathogen population by 0.001 (GPA). Combining this with the above pathogen growth rates gives the following equation for the total population growth of pathogens: dP⁄dt = 1.6 + sgn(PP * 100) + 1/2 - 0.001PA, where sgn(PP * 100) + 1/2 is just a function that returns 1 if PP is greater than 100, and 0 if PP is less than 100.
Two classes of cytokines were used in our model: pro-inflammatory (Pr) and anti-inflammatory (PT). The concentration/population of these cytokines depends directly on PH and on the ratio between the active and suppressive leukocyte populations → Pr = PA + c/PS + c * PH/PmaxM, PT = PS + c/PA + c * PH/PmaxM , where ‘c’ is a smoothing constant used to mitigate the impact of extreme values of PA/PS and ‘M’ is the estimated molarity of the cytokines in the cell-environment.
These cytokine levels also impact PH, PA, and PS: Pr upregulates PH and PA at rate ‘k’, and downregulates PS at the same rate. PT on the other hand downregulates PH and PA at rate ‘k’, and upregulates PS at the same rate. Adding the effects of cytokines on the cell populations gives the following: dH⁄dt = GHPHɸ - DHPH + dr⁄dtk - dT⁄dt
dA⁄dt = (dK⁄dt * dP⁄dt) * Prk
dS⁄dt = (dK⁄dt * dP⁄dt) * PTk
Integration
Our team utilized the Netlogo software, which can be found here. NetLogo is a widely used agent-based modeling software that offers useful tools for developing a simple user-interface for our model. We specifically used version 6.4.0 when coding our tool, but our model is backwards-compatible to at least version 5.0.0 and will be forward-compatible for the foreseeable future.
Using the NetLogo software, the Dry Lab team first set up the underlying cell environment. This was created using “patches,” which are divisions of a two-dimensional space to create a grid-like system. Our model uses a 25 patch x 25 patch environment that wraps vertically and horizontally, meaning that objects moving across one border are transported to the opposite border.
Then, global variables were defined for every growth, death, and differentiation rate mentioned in the above equations. A few software-specific variables were also defined, but those only impact the visualization of our model and have no actual impact on the evolution of the system. One such variable that bears noting is the ‘ticks’ variable. The NetLogo software measures timesteps using ‘ticks,’ which can be defined to represent any fixed amount of ‘real-world’ time. For our model, we defined one tick (timestep) to be 10 minutes. In our results, our time variable is also measured in ticks. It is very important to keep this conversion rate in mind when interpreting the results produced by our model.
The underlying cell layer was defined as a patch variable, which means the agent does not move but can contain many variables indicating its state. To help with the visualization of our model, we added a function that changed the color and intensity of each patch to reflect the relative population of HSPC. The higher the population, the brighter the patch.
Lastly, we added in dynamic plotting, which plots the value of a specific variable at an instantaneous point in time. We also added a way for the user to change viewing modes, allowing them to see the relative concentrations of a particular cell type over time depending on which mode they choose (0,1,2,...7). The images below show snapshots of the cell populations at the same moment in time, and are meant to showcase the different viewing modes:
Using Our Model
All of our work regarding the dry lab model has been documented in this GitHub repository: https://github.com/cfaison/UFlorida2024.git. This is where you can find and download the most current version of our model, as well as view our developmental process.
To use our software, first download the most recent version of our model from the aforementioned GitHub repository. Then open the model using either of the following methods:
- Download and install the NetLogo application to your device: https://ccl.northwestern.edu/netlogo/. Then open our model directly from the application.
- Open the web-based version of NetLogo: https://www.netlogoweb.org/. Then open our model from the browser window.