# Modeling Of Asymmetric Cell Division Biology Essay

**Published:** **Last Edited:**

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Hematopoietic stem cells are distinguished by their differentiation and self-renewal ability. Self-renewal is to maintain the stem cell pool while differentiation on the other hand enables the HSCs to differentiate into later stages of cells. The progenitor cells from the downstream compartments also display both of these abilities. It is yet to be found if the external signal regulates the self-renewal rate or the proliferation rate. Three multi-compartment models which depend on a single feedback mechanism have been proposed. In Model 1, the external signal regulates the proliferation rate, while the self-renewal rate is kept constant. In Model 2, the external signal regulates the self-renewal rate, while the proliferation rate is kept constant. Model 3 is a combination of Model 1 and Model 2 as the external signal regulates both the self-renewal rate and proliferation rate. The main problem lies in determining which of the above models would agree with the system the most. We have employed the ODE 45 solver (4th and 5th order Runge-Kutta formulas) of Matlab to solve these mathematical models. This study highlights that regulation of the fraction of self-renewal rate (Model 2) is the only model that can reach a distinctive absolute positive stable steady state. Furthermore, a smaller amount of successful cell proliferations is essential to reaching this steady state. With the aim of maintaining the stem cell pool in mind, the ratio of self-renewal rate to differentiation has to be higher than the self-renewal rate of downstream compartments and this value has to be more than or equal to 50%. Surprisingly, the equilibrium state of mature cells is independent of the upstream compartments' parameters of dynamics but dependent on parameters of self-renewal rate of hematopoietic stem cells. The model has proven to be valid with respect to leukocyte numbers following hematopoietic stem cell transplantation. This study has explicitly shown that the extrinsic stimulation of the fraction of self-renewal rate is most crucial in the blood cell maturation process.

Introduction

## 1.1 What is systems and computational biology?

Systems biology is a field of study that focuses on the complex interactions in biological systems. It provides a holistic view of how a collection of components interacts differently compared to each component; this is known as emergent properties. Systems biology started as the quantitative modeling of the kinetics of enzyme. It soon expanded to encompass mathematical modeling of population growth and into many other disciplines of studies. Currently there are several main platforms of study; they include phenomics (variation of phenotype changes), genomics (variation in DNA sequences between cells), and transcriptomics (whole cell expression measurements). These platforms identify and quantify the dynamics and interaction within a cell. As there are many possible combinations of results and observations from system biology, a systematic and automated approach must be employed to investigate. This leads to the use of mechanistic modeling that describes biological components into simple or complex series of differential equations. This is where computational biology is used to help speed up the process of simulating mathematical models derived from system biology.

## 1.2 Why do we need it?

The use of computers and the appropriate software enables scientists to simulate the many possible outcomes even before carrying out wet lab experiments. This increases the overall throughput as well as efficiency of research.

## 1.3 What novel insights can it provide, in general?

Through systems and computational biology, we are able to simulate the proliferation and self-renewal of stem cells. This provide a multitude of information as to the process in which cell renew to replace itself as well as differentiate into the mature cells.

## 1.4 What mathematical tools and software are usually adopted in such study?

Mathematical tools such as bifurcation analysis and stability analysis are often adopted in such studies to determine how the choice of parameters would affect the dynamics of the model. Bifurcation analysis is the study of the change of number of steady states when parameters are altered. This is carried out by using an initial guess for certain parameters to simulate the model and eventually comparing the corresponding results with experimental data. The parameters are then adjusted and the whole cycle repeated until the desired results are achieved. After the number of steady states is determined, stability analysis is usually performed. This is done via the Jacobian matrix where the trace and determinant of the matrix are used to verify if a particular steady state is stable or otherwise. In addition, Bode plots and Nyquist plots can be used to further investigate the parameters and their corresponding effects on the stability of the steady states.

The use of Matlab and its huge repository of mathematical tools aid and facilitate the modeling process. It provides the tools to solve ordinary differential equations with initial value conditions. Next to analyze bifurcation and stability, XPPaut provides the mean to evaluate and solve a series of differential equations with the ability to plot and analyze the stability and bifurcation points in the system.

## 1.5 Purpose of choosing the system

Stem cell biology is one of the key potential treatments of degenerative illness, however the mechanism of how HSCs' asymmetric division is controlled/stimulated/regulated has always been an issue that needs to be tackled in stem cell biological arena.

The purpose of this paper is to identify the key factors affecting stem cells therapy, by looking into the underlining mechanism of asymmetric division using computational modeling. By selecting this paper, we hope to give the readers an overview of the possible solutions to this challenge.

## 1.6 Introduction on the biology of the system

Hematopoiesis is a process consisting of multiple steps which generates diverse mature blood cells. A progeny cell, the essential component to hematopoiesis, will either differentiate or replenish the stem cell pool. The multiple roles to differentiate and self-renew is not only limited to HSCs but could also be extended to downstream processes although at a smaller extent. The probability of a daughter cell being identical to its parent cell is equal to self-renewal rate. Conversely, proliferation rate is not equal to self-renewal rate and thus could not be equal to the probability as previously mentioned. Clear evidence shows that asymmetric division is a process that determines the progeny cells' divergent fates during cell division. The research paper demonstrated that slower division rates are closely linked with higher self-renewal rates at a single-cell level. Proof to identify whether the divergent fates of progeny cell are predetermined during cytokinesis or before cytokinesis is still not within reach. Nevertheless, observation by using animal models shows that asymmetric division is valid in the human hematopoietic system. Hypothetically, the fraction of self-renewal rate will be either controlled by external signal or intrinsic mechanism. However, precise clinical data for this question is still unknown.

## 1.7 Detailed description of the biology in our system

Hematopoietic stem cells are found in the bone marrow and representing a quiescent fraction. Through a number of differentiations in 6 compartments, HSCs will eventually transform into diverse mature blood cells.

Figure 1: Process of hematopoiesis.

There is no certain indication of the number of differentiations due to the number of differentiation in each compartment not being fixed. As a result, the number of compartments is considerably lower than the number of differentiations. Cells in upstream compartments replicate at a slower pace (lesser than once per day) as opposed to cells in downstream compartments that undergo faster replication.

After a high-dose of chemotherapy, hematopoiesis will be impaired and leukocyte counts in the peripheral blood will decrease rapidly. As a result, regenerative pressure will rise. Chemotherapy affects every stage of hematopoiesis, hence it is no doubt that the HSC pool will be affected as well. Under normal circumstances, HSC transplantation is performed after high-dose chemotherapy to ensure hematopoietic recovery. Currently, HSC transplantation is no longer a highly experimental procedure, but instead has evolved into a standard therapy for several malignant and hereditary diseases. HSC transplantation can be done through two types of transplantation namely, autologous transplantation (HSC harvested from the patient before chemotherapy) and allogeneic transplantation (HSC harvested from another healthy donor).

Figure 2: Leukocytes per nL versus Days After Transplantation. (The leukocytes counts were obtained from patients with mean Â±SD of all blood composition after transplantation.) High-dose chemotherapy were treated during Days -3 and -2 (Black arrows). Autologous transplantation took place during Day 0 (Grey arrow). Solid line indicate the increment of leukocyte counts after autologous transplantation.

The peripheral blood's leukocyte counts of 20 chemotherapy patients were shown in figure 2. Two days after the high-dose chemotherapy treatment, all patients undergo autologous transplantation. The drop of leukocyte counts from 4 to 10 Ã- 109 leukocytes/L to less than 0.05 Ã- 109 leukocytes/L at Day 7 represents hematopoiesis' retarded effect due to chemotherapy. After which, the leukocyte counts rebound leading to the recovery of the immune system. The leukocyte counts will reach as high as 109 leukocyte/L over an average of 14 days. Nonetheless, the duration for hematopoietic recovery is reliant on progenitor cells which are far more abundant as compared to hematopoietic stem cells.

Hematopoiesis has to be closely regulated by either an external signal or intrinsic mechanism. The mechanism of how HSCs' asymmetric division is controlled / regulated remains elusive in stem cell biology. This issue requires the need of a mathematical modeling paradigm.

The stabilities of the three mathematical models that represent different forms of regulation are evaluated and subsequently compared to clinical data for reproducibility. The 3 models that are termed as Model 1, Model 2, and Model 3 all rely on a single feedback loop.

Figure 3: Model 1

Figure 3 shows the Model 1 as well as the feedback signal that regulates this process. The proliferation rate is regulated directly by the signal while the fraction of self-renewal rate remains constant. When the progenitor cell undergoes asymmetric division, one of the two daughter cells produced would be identical to the parent cell in terms of its properties and is used for replenishment of the stem cell pool. The other daughter cell would differentiate and proceed to the next compartment.

Figure 4: Model 2

Figure 4 shows the Model 2 as well as the feedback signal that regulates this process. The fraction self-renewal rate is regulated directly by the signal while the proliferation rate remains constant. When the progenitor cell undergoes asymmetric division, it is not necessary for one daughter cell to undergo self-renewal while the other differentiates. It has 3 possible outcomes with the first being identical to that of Model 1. The second outcome would be both daughter cells undergoing self-renewal while the last outcome involves both daughter cells proceeding towards differentiation.

Figure 5: Model 3

Figure 5 shows the Model 3 as well as the feedback signal that regulates this process. This model is the combination of Model 1 and Model 2 where the signal regulates both the fraction of self-renewal rate and proliferation rate. Numerical simulation indicates Model 2 is a more effective mechanism than Model 1. Stimulation of asymmetric division via regulation of self-renewal rate provides higher flexibility which in turn would provide faster regeneration as well. Additionally, Model 2 stays in the range that agrees with experimental data on autologous HSC transplantation.

## Materials and Methods

## 2.1 Clinical data for hematopoietic reconstitution

20 patients suffering from advanced multiple myeloma aged 37 to 66 years old, with a median of 44 years old, provided the clinical data involving hematopoietic reconstitution for HSC transplantation (Figure 2). A myeloablative high-dose chemotherapy, with Melphalan 100 mg/m2 body surface, was used on days -3 and -2 to perform the treatment. It was then accompanied by autografting of an unmanipulated HSC transplant, with 2.0 to 10.3 (median of 3.6) CD34+ cells/kg body weight, on day 0. This transplant was harvested prior to the treatment from granulocyte colony-stimulating factor mobilized peripheral blood via leukapheresis. All the above procedures were carried out at the Department of Medicine V in the University of Heidelberg.

## 2.2 Derivation of the multicompartment model

A model was designed to simulate the six compartments, namely long-term repopulating stem cells (with a level defined as c1), short-term repopulating stem cells (c2), multipotent progenitor cells (c3), committed progenitor cells (c4), precursors (c5) and mature cells (c6).

Cell cycle can be defined by ordinary differential equations and the dynamics of hematopoietic cell differentiation and proliferation are described as follows: Cell cycle is assumed to be a well-mixed population, where cells may enter division at a rate p(t), or die at a rate d that can be considered to be a constant factor in time. To simplify the design of the model, the death rate is considered to occur only in the last compartment that consists of mature cells.

For the first compartment containing c1 cells, the flux to mitosis is given by p1(t)c1(t) and the flux of just divided cells from mitosis is given by 2p1(t)c1(t). The daughter cells would either replenish the compartment via self-renewal or differentiate into the next stage of cells and move on to the next compartment. The rate of self-renewal is assumed to be a1(t) where a1(t) has a value ranging from 0 to 1. This means that each daughter cell has the probability of 1-a1(t) of differentiating into the next stage of cells and move on to the next compartment. Since only cell characteristics that are averaged over the entire population are involved in the model, it is required to specify only a certain fraction of the daughter cells, equivalent to 1-a1(t), which would differentiate upon division. Hence, the influx of cells into the cycle of the first stage of differentiation can be defined as 2ai(t)pi(t)ci(t) and the resultant rate of change of c1(t) can be determined from the following equation:

Now consider the subsequent stages and their respective compartments i, where i = 2-5, with ci cells. The differentiation of the cells can be described as a flow of cells from the previous compartment i-1 to the current compartment i. Cell replication in ci can be taken to be similar as that of the first compartment c1, meaning an influx of (2ai(t)-1)pi(t)ci(t) cells per unit time. In addition, there are progeny cells transported from the previous compartment ci-1 at a rate of 1-ai-1(t). Hence, there is a further increase in the total cell number in compartment i and it is given by:

For the sixth compartment, the mature cells do not proliferate anymore and its increase is only provided by the flux of differentiated cells from the previous compartment, namely compartment five. The cell number is reduced only by the death of the mature cells dc6(t) and it can be represented by the following formula:

The proliferation and asymmetric cell divisions rates, defined as pi(t) and ai(t) respectively, can be time-dependent functions. In the following section of the paper, we would take them as either being constants or functions of a feedback signal that can be generated by a mechanism that is described by a separate differential equation.

## 2.3 Model of the feedback signal

Cytokines are extracellular signaling molecules that aid in the control of the dynamics of cell differentiation and proliferation. The exact nature of this process is still not fully understood. However, we can assume that cytokines are secreted by specialized cells at a maximum rate of Î± and this secretion is regulated by a negative feedback mechanism that is determined by the concentration of the secreted cytokines, ÂµS. This means that the production of cytokines is in fact a self-limiting process. In addition, the supply of cytokines is further regulated by mechanisms that are sensitive to the number of mature cells present, Î²Sc6. This means that as the number of mature cells increases, the production of the signaling factor would be decreased. From these assumptions, we can construct the following differential equation for the dynamics of the cytokines:

Substituting the parameters s = Âµ/Î± âˆ™ S and k = Î²/Âµ into the above equation, we would obtain the resultant equation ds/dt = Âµ(1-s-ksc6), which allows us to understand how the relative intensity of the feedback signal varies with time. However, since the production of cytokines is very fast as compared to the rates of cell differentiation and proliferation, where it may take up to days, the level of cytokines would approach a steady state very quickly. Hence, by applying a quasi-steady state approximation, we would be able to obtain the following equation that allows us to determine the feedback signal intensity:

This resultant model produces the assumption that the feedback signal is at a maximum when mature cells are completely absent from the system. On the other hand, as the number of mature cells increases, the signal would decrease asymptotically to zero. Also, k is taken to have the value of 1.6 x 10-10 in the subsequent numerical simulations.

## 2.4 Description of mathematical equations

## Model 1

For the first model, the feedback signal from the cytokines is used to regulate the proliferation rates (pi) of all the six compartments while the asymmetry of cell division is left unchanged. All of the hypotheses mentioned above are used to construct a system of six ordinary differential equations for the six respective compartments (c1, c2, ..., c6) with each of them being a function of time:

## Model 2

For the second model, the feedback signal from the cytokines is used to regulate the asymmetry of cell division, which is the ratio of the rate of self-renewal to the rate of differentiation, ai. However, unlike Model 1, the proliferation rate is left unchanged:

## Model 3

For the third model, the feedback signal from the cytokines is used to regulate both the parameters, namely the proliferation rate and the asymmetry of cell division, at the same time:

## 2.5 Overview

## Model

## Rate of Self Renewal

## Proliferation Rate

## 1

Constant

Signal Regulated

## 2

Signal Regulated

Constant

## 3

Signal Regulated

Signal Regulated

To test the three different hypotheses related to the mechanisms of regulating the ratio of rate of differentiation by feedback signaling factors to rate of self renewal, the team has designed three sets of differential equations. These three sets of mathematical modeling takes into consideration the reconstitution of leukocytes after hemopoietic stem cell transplant treatment after chemotherapy. The initial conditions used for the 6 compartments are as follows, C1(0) = 105, C2(0) = 106, C3(0) = 107, C4(0) = 107 C5(0) = 107, C6(0) = 107. Using the initial conditions, the numbers of matured cells were plotted as a function of time (0 to 100 days).

## Results

## Model 1

This model describes the case where proliferation rate is signal regulated, where the increase in signaling molecules affects the proliferation rate, while the rate of self renewal is kept constant. In this model, the rate of replication is observed to increase with the decreasing number of mature cells. This leads to the stabilization of cell population at an equilibrium point. The equation for compartment 1 used in this model indicates that the system has a non trivial equilibrium when 50% of the long term repopulating stem cells replenishes the short term repopulating stem cells, this leads to non isolated equilibrium points. The steady state assumption is dependent on the initial population size of HSC.

## Model 2

This model describes the hypothesis that during aplasia, signaling molecules are secreted which influences the ratio of differentiation rate versus self renewal rate. High levels of signaling molecules increases and speed up self renewal, while low levels of signaling molecules encourage differentiation, forcing the cells to differentiate to the next compartment. In this model, a positive steady state exist when the hematopoietic stem cell self renewal (a1Â) is larger than the other component's self renewal rates such as (a2 a3 a4 a5 a6). Therefore it is critical that the steam cell progeny pool has the largest population as it contributes downstream into the other compartments. In an event where any compartment other then the HSC population has a significantly larger amount, it will take over the stem cell functions causing all upstream compartments including the HSC compartment to extinct with time, this observation is consistent with other authors. Interestingly, the steady state level of mature cell population depends on the self renewal rate (a1Â) of the HSC population as seen from the equation used for compartment 6, . From linear stability analysis, the equilibrium exists in a local asymptotical manner, where as other steady states such as trivial and semi-trivial are not stable. Using numerical simulation, the steady state obtain is globally stable and converges. In comparison, the steady state for model 1 does not depend on initial cell numbers, and therefore has the ability to equalize any lost of cell numbers. This model has an exceptionally close response to that observed in clinical data.

## Model 3

In this model, the feedback signals regulate both the proliferation rate and rate of self renewal. The behavior of this model is qualitatively similar to that of model 2. This model has a unique steady state which is asymptotically stable when tested using numerical simulations as well as linear stability analysis.

## 3.1 Matlab Codes for Modeling

To solve the differential equations mentioned above, we would first have to process them into string of codes which enables us to solve them using matlab. The following codes were modeled according to the three different models mentioned in the paper with the appropriate parameters extracted from the paper. Comments have been included in the codes to describe each section and explain the purpose of that section.

## Code for Model 1

functiondcdt = Model_1(t, c)

%parameters documented in the paper

pa=0.173; %p1

pb=0.231; %p2

pc=0.347; %p3

pd=0.693; %p4

pe=1.386; %p5

aa=0.5; %a1

ab=0.475; %a2

ac=0.475; %a3

ad=0.475; %a4

ae=0.41; %a5

k=1.6*10.^-10; %constant k

f=0.3; %d

dcdt = zeros (6,1);

%following are the 6 compartment ODE

dcdt(1)= (2*aa-1)*pa*(1/(1+k*c(6)))*c(1);

dcdt(2)= ((2*ab-1)*pb*(1/(1+k*c(6)))*c(2))+(2*(1-aa)*pa*(1/(1+k*c(6)))*c(1));

dcdt(3)= ((2*ac-1)*pc*(1/(1+k*c(6)))*c(3))+(2*(1-ab)*pb*(1/(1+k*c(6)))*c(2));

dcdt(4)= ((2*ad-1)*pd*(1/(1+k*c(6)))*c(4))+(2*(1-ac)*pc*(1/(1+k*c(6)))*c(3));

dcdt(5)= ((2*ae-1)*pe*(1/(1+k*c(6)))*c(5))+(2*(1-ad)*pd*(1/(1+k*c(6)))*c(4));

dcdt(6)= (2*(1-ae)*pe*(1/(1+k*c(6)))*c(5))-(f*c(6));

## Code for Model 2

functiondcdt = Model_2(t, c)

%parameters documented in the paper

pa=0.173/2; %p1

pb=0.231/2; %p2

pc=0.347/2; %p3

pd=0.693/2; %p4

pe=1.386/2; %p5

aa=0.7; %a1

ab=0.65; %a2

ac=0.65; %a3

ad=0.65; %a4

ae=0.55; %a5

k=1.6*10.^-10; %constant k

f=0.3; %d

dcdt = zeros (6,1);

%following are the 6 compartment ODE

dcdt(1)= (2*aa*(1/(1+k*c(6)))-1)*pa*c(1);

dcdt(2)= ((2*ab*(1/(1+k*c(6)))-1)*pb*c(2))+(2*(1-aa*(1/(1+k*c(6))))*pa*c(1));

dcdt(3)= ((2*ac*(1/(1+k*c(6)))-1)*pc*c(3))+(2*(1-ab*(1/(1+k*c(6))))*pb*c(2));

dcdt(4)= ((2*ad*(1/(1+k*c(6)))-1)*pd*c(4))+(2*(1-ac*(1/(1+k*c(6))))*pc*c(3));

dcdt(5)= ((2*ae*(1/(1+k*c(6)))-1)*pe*c(5))+(2*(1-ad*(1/(1+k*c(6))))*pd*c(4));

dcdt(6)= (2*(1-ae*(1/(1+k*c(6))))*pe*c(5))-(f*c(6));

## Code for Model 3

functiondcdt = Model_3(t, c)

%parameters documented in the paper

pa=0.173; %p1

pb=0.231; %p2

pc=0.347; %p3

pd=0.693; %p4

pe=1.386; %p5

aa=0.7; %a1

ab=0.65; %a2

ac=0.65; %a3

ad=0.65; %a4

ae=0.55; %a5

k=1.6*10.^-10; %constant k

f=0.3; %d

dcdt = zeros (6,1);

%following are the 6 compartment ODE

dcdt(1)= (2*aa*(1/(1+k*c(6)))-1)*pa*c(1)*(1/(1+k*c(6)));

dcdt(2)= ((2*ab*(1/(1+k*c(6)))-1)*pb*c(2)*(1/(1+k*c(6))))+(2*(1-aa*(1/(1+k*c(6))))*pa*c(1)*(1/(1+k*c(6))));

dcdt(3)= ((2*ac*(1/(1+k*c(6)))-1)*pc*c(3)*(1/(1+k*c(6))))+(2*(1-ab*(1/(1+k*c(6))))*pb*c(2)*(1/(1+k*c(6))));

dcdt(4)= ((2*ad*(1/(1+k*c(6)))-1)*pd*c(4)*(1/(1+k*c(6))))+(2*(1-ac*(1/(1+k*c(6))))*pc*c(3)*(1/(1+k*c(6))));

dcdt(5)= ((2*ae*(1/(1+k*c(6)))-1)*pe*c(5)*(1/(1+k*c(6))))+(2*(1-ad*(1/(1+k*c(6))))*pd*c(4)*(1/(1+k*c(6))));

dcdt(6)= (2*(1-ae*(1/(1+k*c(6))))*pe*c(5)*(1/(1+k*c(6))))-(f*c(6));

In order to display the values of each compartment in graphical representation, and compare it across the three different models, the following code was used to solve and plot the graphs. We used the ODE 45 solver in matlab as it is an explicit Runge-Kutta 4th / 5th order method and is efficient in our application. We included the initial conditions mentioned in the paper, as well as customized the plot function to provide us with the ideal graphs.

## Code for Compartment 1

## Code for Compartment 2

functionCompartment_1

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,1),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,1),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,1),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

functionCompartment_2

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,2),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,2),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,2),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

## Code for Compartment 3

## Code for Compartment 4

functionCompartment_3

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,3),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,3),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,3),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

functionCompartment_4

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,4),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,4),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,4),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

## Code for Compartment 5

## Code for Compartment 6

functionCompartment_5

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,5),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,5),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,5),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

functionCompartment_6

%using ode45 solver with parameters documented in the paper

%time axis from 0 to 100 days

%initial cell population

[t,c] = ode45('Model_1',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,6),'b');

Hold all;

[t,c] = ode45('Model_2',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,6),'r');

Hold all;

[t,c] = ode45('Model_3',[0 100],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c(:,6),'black');

%plotting parameters to customize graph

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Model 1', 'Model 2', 'Model 3');

Using the above matlab codes, we were able to generate the following graphs. The graphs obtained were very similar in terms of its shape and the rate in which each compartment takes to reach steady state. Interesting, the steady state of compartments 1, 2 and 3 generated from our code do not yield similar maximum cell numbers compared to the paper, we will discuss this observation in the summary portion of this report. However, we were able to achieve the same graphs for the rest of the figures.

Figure 6: Modeling of hematopoietic reconstitution after HSC transplant treatment.

(A-F; Model 1: Blue line; Model 2: Red line; Model 3: Black line.)

## 3.2 Additional results (Not in research paper)

Briefly mentioned in the paper was the idea of a model that included the outflux of hematopoietic stem cells. This outflux could be due to the death of HSC or the transition of state into quiescence, we denote it as (d1) and the following equation is obtained.

Although the paper did not provide the results of this modeling, our team managed to solve it and came out with the results. When we solve the above equation, we obtain the graph shown in figure 7. An overall decrease in the maximum cell population is observed when compared to model 1. In addition, massive decline in cell population is observed at the 600 days mark. The trajectory calculated in this alternate model 1 does not match the hematopoietic reconstitution curve as the growth of the matured cells is too slow, therefore this model is rejected.

## Code for Alternate Model 1

functiondcdt = Model_1_alternative(t, c)

%parameters documented in the paper

pa=0.173; %p1

pb=0.231; %p2

pc=0.347; %p3

pd=0.693; %p4

pe=1.386; %p5

aa=0.5; %a1

ab=0.475; %a2

ac=0.475; %a3

ad=0.475; %a4

ae=0.41; %a5

k=1.6*10.^-10; %constant k

f=0.3; %d

dcdt = zeros (6,1);

%following are the 6 compartment ODE

dcdt(1)= (2*aa*-1)*pa*(1/(1+k*c(6)))*c(1)-(f*c(1));

dcdt(2)= ((2*ab-1)*pb*(1/(1+k*c(6)))*c(2))+(2*(1-aa)*pa*(1/(1+k*c(6)))*c(1));

dcdt(3)= ((2*ac-1)*pc*(1/(1+k*c(6)))*c(3))+(2*(1-ab)*pb*(1/(1+k*c(6)))*c(2));

dcdt(4)= ((2*ad-1)*pd*(1/(1+k*c(6)))*c(4))+(2*(1-ac)*pc*(1/(1+k*c(6)))*c(3));

dcdt(5)= ((2*ae-1)*pe*(1/(1+k*c(6)))*c(5))+(2*(1-ad)*pd*(1/(1+k*c(6)))*c(4));

dcdt(6)= (2*(1-ae)*pe*(1/(1+k*c(6)))*c(5))-(f*c(6));

## Code to plot alternative model

functionModel_1_Alt_all_Compartment

%plot of the alternate model 1 graph

[t,c] = ode45('Model_1_alternative',[0 1000],[100000, 1000000, 10000000, 10000000, 10000000, 10000000]);

plot(t,c);

ylabel('Number of Cells');

set(gca,'YMinorTick','on');

xlabel('Days');

Legend('Compartment 1', 'Compartment 2', 'Compartment 3','Compartment 4','Compartment 5','Compartment 6');

Figure 7: Alternative modeling of hematopoietic reconstitution when outflux is considered (left), compared to the original model1.

Overall, our team was able to repeat all of the results mentioned in the paper using the tools provided in Matlab and at the same time produced additional results that were not covered in the paper.

## Summary

## 4.1 Significance of simulation

In this modern era, where science and medicine is advanced, incurable diseases such as acute leukemia and hereditary immunodeficiency still exist and have been taking a toll on mankind. These diseases are constantly being researched on but to no avail thanks to our complex biological system. And thus the need to understand this biological system arises. Stem cell therapies have been the new "in-thing" in the medical arena. By breaking frontiers in various treatments of animals, it has been considered as the pioneer of a new breed of treatments. According to Baillie (2008), adult stem cells were harvested from a dog's fat cells and were used to treat its arthritic condition. Such successful treatments involving animals call for a greater understanding of stem cell therapies in order to apply them to humans. Hematopoietic stem cell (HSC) being the very essence of these therapies play an important role in the process of hematopoiesis. Being a key player in this process, it is essential to understand the HSC and this is what this report has been focusing on. Through pointing the direction at the extrinsic regulation of the self-renewal rate of HSC which serves as the key process of hematopoiesis, this report aims to provide the reader an insight into achieving efficient repopulation. Simply put, the significance of the simulation in this report is to give an overview into the process of hematopoiesis and how different variables could affect the outcome of this process. It is believed that through this understanding efficient and effective treatments could be developed for humans.

## 4.2 Difficulties encountered

A great mind once said, "To be successful you must accept all challenges that come your way. You can't just accept the ones you like". With reference to the quote, there were indeed some challenges that we came across and attempted to solve while simulating the model. As previously addressed, matlab software was employed to simulate the model. Though the general trends of the graphs were achieved, the physical quantities did not match i.e. cell numbers of our model did not match with the cell numbers of the actual model proposed by the authors in the research paper. Considering that it might be an error as a result of choosing the wrong ode solver, we simulated the model with other ode solvers i.e. ode 23, ode 113 and so on. It soon came to our attention that the results further deviated from that in the research paper with varied ode solvers and did not make sense. Subsequently, we stuck to ode 45 as that was the best possible option to solve the set of differential equations. This meant the problem had to lie elsewhere. Upon further exploration, we realized that the problem was with the initial conditions that we based on to solve the differential equations. Initial conditions were only specified for the first three compartments by the research paper and that we assumed certain values for the subsequent compartments. We attempted at playing around with the numbers but to no avail. The exact values of the model could not be obtained. On logical thinking, we realized that the trends of the graph were more important as compared to the exact values which led to our main objective in the determining the mechanism of regulation that best suits hematopoiesis. Numbers did not mean anything to a certain extent as priority was to understand the system itself.

## 4.3 Future considerations (Not in Research Paper)

Long story short, this model focused on two variables namely, proliferation rate and self-renewal rate and how they effected changes to hematopoiesis. By taking a step further, future models could focus on aspects that may affect self-renewal rate and proliferation rate. During our desk research, it came to light that multiple factors such as tumor suppression genes, transcription factors and cell polarity could influence self-renewal or proliferation. Dysregulation of these factors could spell disaster in the form of cancer and developmental disorders. In addition, these factors could provide additional parameters that will demonstrate how defects in self-renewal and proliferation could affect the process of asymmetry of cell divisions. Through the extension of the current model, future models could aid in identifying stem cells and their environment and thus will lead to a new perspective of control of self-renewal and proliferation.

This model is based upon autologous transplants and has incorporated data from autologous transplantations conducted in mice. For starters, autologous transplant refers to transplanted stem cells that come from your own blood or bone marrow. The shortcoming of this model is that it could not be applied to allogenic transplants where donor tissue does not originate from the patient itself considering that other external variables will be present to affect hematopoiesis. In an example pertaining to blood stem cell transplant, only five out of six HLA antigens were used as perfect 6 matches were hard to achieve. This ultimately led to issues such as graft rejection and graft versus host disease. Hence, future models could be developed to incorporate allogenic transplantation which is the most feasible mode of transplantation and along the modeling process issues such as mismatched tissue could be accounted for. These new models could pave the path in understanding how tissues or cells of different origins could affect hematopoiesis.