Analyzing the birth-death model of Oncostreams in Glioma, and the effects of Cytochalasin D treatment

In a birth-death processes the appearance of new individual bacterial cells, as well as the disappearance of existing individual bacterial cells play a very important role. Looking into simple birth and simple death processes will allow us to build a basis to expand for further applications going forward. We introduce important values and notation as follows below.

	$\displaystyle\lambda$	$\displaystyle\rightarrow\text{birth rate}$
	$\displaystyle\mu$	$\displaystyle\rightarrow\text{death rate}$
	$\displaystyle N(t)\text{ and }n$	$\displaystyle\rightarrow\text{population}$
	$\displaystyle\lambda N(t)\Delta t$	$\displaystyle\rightarrow\text{chance of +1 (1 birth)}$
	$\displaystyle\mu N(t)\Delta t$	$\displaystyle\rightarrow\text{chance of -1 (1 death)}$
	$\displaystyle 1-(\lambda+\mu)N(t)\Delta t$	$\displaystyle\rightarrow\text{chance of 0 (no change)}$
	$\displaystyle a$	$\displaystyle\rightarrow\text{Initial Population of Glioma Cells}$

Given a birth rate of $\lambda$ of a population of size $N(t)$ with $n$ individuals, we can now create a Kolmogorov equation, which can be solved and used to model this phenomena.

(1)

$\displaystyle p_{n}(t+\Delta t)$

$\displaystyle=p_{n-1}(t)\lambda(n-1)\Delta t+\mu(n+1)p_{n+1}(t)\Delta t-\left(\lambda+\mu\right)n\Delta tp_{n}(t)$

Gives us the difference equation:

(2)

$\displaystyle\frac{dp_{n}(t)}{dt}$

$\displaystyle=\lambda(n-1)p_{n-1}(t)+\mu(n+1)p_{n+1}(t)-(\lambda+\mu)np_{n}(t)$

Notice that if the population becomes $0$ then no births can occur, so we are only concerned about the process starting with a non-zero population suppose it’s $N(0)=a$, leading to an initial condition:

$$p_{a}(0)=1$$

Some things to note here:

	$\displaystyle\mu$	$\displaystyle=0:$
	$\displaystyle\frac{dp_{n}(t)}{dt}$	$\displaystyle=\lambda(n-1)p_{n-1}(t)-\lambda np_{n}(t)$
	$\displaystyle\lambda$	$\displaystyle=0:$
	$\displaystyle\frac{dp_{n}(t)}{dt}$	$\displaystyle=\mu(n+1)p_{n+1}(t)-\mu np_{n}(t)$

Now using steps shown in 6

we get the probability generating function for $\lambda=\mu$:

(3)

$\displaystyle G(s,t)$

$\displaystyle=\left(\frac{1-(\lambda t-1)(s-1)}{1-\lambda t(s-1)}\right)^{a}$

and for the probability generating function for $\lambda\neq\mu$

(4)

$\displaystyle G(s,t)=\left(\frac{\mu\omega(s,t)-1}{\lambda\omega(s,t)-1}\right)^{a}\text{ , where }\omega(s,t)=\frac{(s-1)e^{(\lambda-\mu)t}}{\lambda s-\mu}$

And we can then represent the value of $p_{n}(t)$ as a specific distribution:

	$\displaystyle p_{n}(t)$	$\displaystyle=\frac{(\lambda t)^{n-1}}{(1+\lambda t)^{n+1}}\text{, }n\geq 1$
	$\displaystyle p_{0}(t)$	$\displaystyle=\frac{\lambda t}{1+\lambda t}$

We can follow a similar derivation process for $\lambda\neq\mu$ for $a=1$, which results in:

	$\displaystyle p_{n}(t)$	$\displaystyle=(1-\alpha)(1-\beta)\beta^{n-1}\text{, }\alpha=\frac{\mu(e^{(\lambda-\mu)t}-1)}{\lambda e^{(\lambda-\mu)t}-\mu}\text{, }\beta=\frac{\lambda(e^{(\lambda-\mu)t}-1)}{\lambda e^{(\lambda-\mu)t}-\mu}\text{, }n\geq 1$
	$\displaystyle p_{0}(t)$	$\displaystyle=\alpha$

We now introduce a Population model of the Cytochalasin D treatment and also a model for the concentration of the Cytochalasin D used to help dismantle the oncostreams to slow the glioma invasiveness. And analyze how and whether the growth of glioma changes.

Drug Background information: Cytochalasin D is a drug that is known to impede actin polymerization, and was observed to disrupt the structural organization of oncostreams in GFP+ NPA glioma cells. According to [HUANG] the half life of natural Cytochalasin D is 10 minutes and liposomal Cytochalasin D is 4 hours (240 minutes). Because the paper doesn’t mention if it uses natural or liposomal Cytochalasin D, we will analyze our model using both types.

	$\displaystyle\gamma$	$\displaystyle\rightarrow\text{Decay rate of Cytochalasin D}$
	$\displaystyle\beta$	$\displaystyle\rightarrow\text{Constant used to scale the effect of the treatment concentration}$
	$\displaystyle C(t)$	$\displaystyle\rightarrow\text{Concentration of the Cytochalasin D drug administered}$
	$\displaystyle(\mu+C(t))N(t)\Delta t$	$\displaystyle\rightarrow\text{$-1$ cell implying the death of a cell}$
	$\displaystyle 1-(\lambda+\mu+C(t))N(t)\Delta t$	$\displaystyle\rightarrow\text{No change in the net population or growth of the cells}$

For the simplicity of the derivation process, we can assign $\beta=1$. But later, we will change it to be $\beta=\frac{1}{100}$, for reasons explained in 3.1. This will not cause any problems with our derivation because $\beta$ is a constant.

We can now introduce elementary models to model the Population and Concentration of the Cytochalasin D drug. We are assuming that the treatment is being administered at a constant rate in our paper:

(5)		$\displaystyle\frac{dC(t)}{dt}$	$\displaystyle=-\gamma C(t)$
(6)		$\displaystyle\frac{dN(t)}{dt}$	$\displaystyle=\lambda N(t)-\mu N(t)-C(t)N(t)$

We now solve the ODEs to get a more accurate form of the $N(t)$ and $C(t)$ equations:

	$\displaystyle\frac{dC(t)}{dt}$	$\displaystyle=-\gamma C(t)$
	$\displaystyle\int\frac{dC(t)}{C(t)}$	$\displaystyle=\int-\gamma dt$
	$\displaystyle\ln{C(t)}$	$\displaystyle=-\gamma t+A$
	$\displaystyle\text{Using the initial condition of }C(0)$	$\displaystyle=C_{0}:$
	$\displaystyle C(t)$	$\displaystyle=C_{0}e^{-\gamma t}$

Using the this result we can solve for the population function. We can use the assumption that $N(0)=a$, where $a$ is the initial population of glioma cells.

	$\displaystyle\frac{dN(t)}{dt}$	$\displaystyle=\lambda N(t)-\mu N(t)-C_{0}e^{-\gamma t}N(t)$
	$\displaystyle\frac{dN(t)}{dt}$	$\displaystyle=N(t)\left[\lambda-\mu-C_{0}e^{-\gamma t}\right]$
	$\displaystyle\int\frac{dN(t)}{N(t)}$	$\displaystyle=\int[\lambda-\mu-C_{0}e^{-\gamma t}]dt$
	$\displaystyle\ln{N(t)}$	$\displaystyle=\lambda t-\mu t+\frac{C_{0}}{\gamma}e^{-\gamma t}+B$
	$\displaystyle N(t)$	$\displaystyle=\frac{a}{e^{\frac{C_{0}}{\gamma}}}e^{t(\lambda-\mu)+\frac{C_{0}}{\gamma}e^{-\gamma t}}$

Using these values and derivations, we can derive a Kolmogorov Equation:

(7)

$$\frac{dp_{n}(t)}{dt}=\lambda(n-1)p_{n-1}(t)+[\mu+C_{0}e^{-\gamma t}](n+1)p_{n+1}(t)-[(\lambda+\mu)+C_{0}e^{-\gamma t}]np_{n}(t)$$

Assume that we are given the initial condition $p_{a}(0)=1$ with only one singular cell at the beginning of time in the cell population. Note that you can have different starting cell values, however it is classic to start with a singular cell body.

This is an example of the Quasi-Steady Stationary State Approximation, because for small $\gamma$, over short intervals, $C(t)$ changes slowly allowing us to treat the concentration of the treatment cells as approximately linear over these intervals.

Some notes about the terms we have appearing in the Kolmogorov Equation above:

• •

  Birth rate $\lambda$: The term $\lambda(n-1)p_{n-1}(t)$ represents the rate at which the population increases by one unit (birth). The coefficient $\lambda(n-1)$ indicates that the rate is proportional to the current population size minus one.

• •

  Death rate $\mu$: The term $\mu(n+1)p_{n+1}(t)$ represents the rate at which the population decreases by one unit (death). The coefficient $\mu(n+1)$ indicates that the rate is proportional to the current population size plus one.

• •

  Treatment term $C_{0}e^{-\gamma t}(n+1)$: This term represents a treatment or decay process where the rate is influenced by a treatment concentration that decays exponentially over time. $C_{0}$ is the initial concentration and $\gamma$ is the decay rate of the treatment.

• •

  Balance Term: The term $-[(\lambda+\mu)+C_{0}e^{-\gamma t}]np_{n}(t)$ ensures the probability is conserved and accounts for the transitions out of the state $n$.

We can now start solving for the probability generating function, which we will define as $G(s,t)$.

$\displaystyle\begin{split}\frac{dp_{n}(t)}{dt}&=\lambda(n-1)p_{n-1}(t)+[\mu+C_{0}e^{-\gamma t}](n+1)p_{n+1}(t)-[(\lambda+\mu)+C_{0}e^{-\gamma t}]np_{n}(t)\\ \sum_{n=0}^{\infty}s^{n}\frac{dp_{n}(t)}{dt}&=\sum s^{n}[\lambda(n-1)p_{n-1}(t)+[\mu+C_{0}e^{-\gamma t}](n+1)p_{n+1}(t)-[(\lambda+\mu)+C_{0}e^{-\gamma t}]np_{n}(t)]\\ \implies&\sum s^{n}(\lambda(n-1)p_{n-1}(t))+\sum s^{n}(\mu(n+1)p_{n+1}(t))+\sum s^{n}(C_{0}e^{-\gamma t}(n+1)p_{n+1}(t))-\\ &-\sum s^{n}(\lambda+\mu)np_{n}(t)-\sum s^{n}(C_{0}e^{-\gamma t})np_{n}(t)\\ \implies&\lambda s^{2}\sum s^{n-2}(n-1)p_{n-1}(t)+\mu\sum s^{n}(n+1)p_{n+1}(t)+C_{0}e^{-\gamma t}\sum s^{n}(n+1)p_{n+1}(t)-\\ &-(\lambda+\mu)s\sum s^{n-1}np_{n}(t)-C_{0}e^{-\gamma t}s\sum s^{n-1}np_{n}(t)\\ \implies&\lambda s^{2}\frac{d}{ds}\left[\sum s^{n-1}p_{n-1}(t)\right]+\mu\frac{d}{ds}\left[\sum s^{n+1}p_{n+1}(t)\right]+C_{0}e^{-\gamma t}\frac{d}{ds}\left[\sum s^{n+1}p_{n+1}(t)\right]-\\ &-(\lambda+\mu)s\frac{d}{ds}\left[\sum s^{n}p_{n}(t)\right]-C_{0}e^{-\gamma t}s\frac{d}{ds}\left[\sum s^{n}p_{n}(t)\right]\\ \implies\frac{\partial G(s,t)}{\partial t}&=(\lambda s^{2}+\mu+C_{0}e^{-\gamma t}-(\lambda+\mu)s-C_{0}e^{-\gamma t}s)\frac{\partial G(s,t)}{\partial s}\\ \frac{\partial G(s,t)}{\partial t}&=(\lambda s^{2}-s(\lambda+\mu+C_{0}e^{-\gamma t})+(\mu+C_{0}e^{-\gamma t}))\frac{\partial G(s,t)}{\partial s}\end{split}$

As this is a first order partial differential equation, we can solve it easily using Mathematica (code in Appendix 7 ). The final result is as follows:

(8)

$$G(s,t)=C_{1}\left[t-\frac{e^{\gamma t}\log{\left[C_{0}(s-1)+e^{\gamma t}(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}{C_{0}+e^{\gamma t}(\lambda+\mu)}\right]$$

Where $C_{1}$ is an integration constant that can be solved for using the initial condition $G(s,0)=s^{a}$.

To find $C_{1}$, we use the initial condition $G(s,0)=s^{a}$. Plugging $t=0$ into the expression for $G(s,t)$:

$$G(s,0)=C_{1}\left[0-\frac{e^{\gamma\cdot 0}\log{\left[C_{0}(s-1)+e^{\gamma\cdot 0}(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}{C_{0}+e^{\gamma\cdot 0}(\lambda+\mu)}\right]$$

Simplifying this expression:

$$G(s,0)=C_{1}\left[-\frac{\log{\left[C_{0}(s-1)+(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}{C_{0}+\lambda+\mu}\right]$$

Given $G(s,0)=s^{a}$, we can set this equal to the simplified expression above:

$$s^{a}=C_{1}\left[-\frac{\log{\left[C_{0}(s-1)+(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}{C_{0}+\lambda+\mu}\right]$$

To isolate $C_{1}$, we rearrange the equation:

$$C_{1}=\frac{s^{a}(C_{0}+\lambda+\mu)}{-\log{\left[C_{0}(s-1)+(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}$$

Thus, the probability generating function $G(s,t)$ can be written as:

$$G(s,t)=\frac{s^{a}(C_{0}+\lambda+\mu)}{-\log{\left[C_{0}(s-1)+(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}\left[t-\frac{e^{\gamma t}\log{\left[C_{0}(s-1)+e^{\gamma t}(-\lambda^{2}-\mu+s(\lambda+\mu))\right]}}{C_{0}+e^{\gamma t}(\lambda+\mu)}\right]$$

Once we have the probability generating function $G(s,t)$, the next step is to use it to derive the moments or the probabilities $p_{n}(t)$. Here’s how we can proceed:

The first moment, or the mean of the distribution $N(t)$, can be found by differentiating $G(s,t)$ with respect to $s$ and then setting $s=1$.

$$\mathbb{E}[X(t)]=\left.\frac{\partial G(s,t)}{\partial s}\right|_{s=1}$$

Where $X(t)$ is treated as a random variable. Using Mathematica, code in the Appendix 8 we see that the mean value is:

$$-\frac{a(C_{0}+\lambda+\mu)\left(t-\frac{e^{\gamma t}\log\left(\left(\lambda-\lambda^{2}\right)e^{\gamma t}\right)}{C_{0}+(\lambda+\mu)e^{\gamma t}}\right)}{\log\left(\lambda-\lambda^{2}\right)}+\frac{C_{0}+\lambda+\mu}{\left(\lambda-\lambda^{2}\right)\log\left(\lambda-\lambda^{2}\right)}+\frac{(C_{0}+\lambda+\mu)^{2}\left(t-\frac{e^{\gamma t}\log\left(\left(\lambda-\lambda^{2}\right)e^{\gamma t}\right)}{C_{0}+(\lambda+\mu)e^{\gamma t}}\right)}{\left(\lambda-\lambda^{2}\right)\log^{2}\left(\lambda-\lambda^{2}\right)}$$

To compute the second moment, we need to evaluate the second derivative of the probability generating function $G(s,t)$ with respect to $s$ and then set $s=1$.

$$\textit{Var}[X(t)^{2}]=\left.\frac{\partial^{2}G(s,t)}{\partial s^{2}}\right|_{s=1}$$

Let’s compute the second derivative of $G(s,t)$ using the Mathematica code highlighted in 9:

$$-2a(C_{0}+\lambda+\mu)\left(-\frac{(C_{0}+\lambda+\mu)\left(t-\frac{e^{\gamma t}\log\left(\left(\lambda-\lambda^{2}\right)e^{\gamma t}\right)}{C_{0}+(\lambda+\mu)e^{\gamma t}}\right)}{\left(\lambda-\lambda^{2}\right)\log^{2}\left(\lambda-\lambda^{2}\right)}-\frac{1}{\left(\lambda-\lambda^{2}\right)\log\left(\lambda-\lambda^{2}\right)}\right)-\\ -\frac{(a-1)a(C_{0}+\lambda+\mu)\left(t-\frac{e^{\gamma t}\log\left(\left(\lambda-\lambda^{2}\right)e^{\gamma t}\right)}{C_{0}+(\lambda+\mu)e^{\gamma t}}\right)}{\log\left(\lambda-\lambda^{2}\right)}-\\ -(C_{0}+\lambda+\mu)\Bigg{(}\frac{2(C_{0}+\lambda+\mu)}{\left(\lambda-\lambda^{2}\right)^{2}\log^{2}\left(\lambda-\lambda^{2}\right)}+\left(\frac{2(C_{0}+\lambda+\mu)^{2}}{\left(\lambda-\lambda^{2}\right)^{2}\log^{3}\left(\lambda-\lambda^{2}\right)}+\frac{(C_{0}+\lambda+\mu)^{2}}{\left(\lambda-\lambda^{2}\right)^{2}\log^{2}\left(\lambda-\lambda^{2}\right)}\right)\cdot\\ \cdot\left(t-\frac{e^{\gamma t}\log\left(\left(\lambda-\lambda^{2}\right)e^{\gamma t}\right)}{C_{0}+(\lambda+\mu)e^{\gamma t}}\right)+\frac{e^{\gamma(-t)}\left(C_{0}+(\lambda+\mu)e^{\gamma t}\right)}{\left(\lambda-\lambda^{2}\right)^{2}\log\left(\lambda-\lambda^{2}\right)}\Bigg{)}$$

One thing to note here, finding a correlation between the First and Second moment proved to be too hectic and no analysis relating the two has been done. Further work on this can be done in related or future work.

Due to the complexity of finding a pattern and multiple derivatives that are needed in finding the probability formula $p_{n}(t)$, we pursue a more numerical approach and numerically solve the problem. The code is included in the Appendix 11, and the motivation for this code was derived from here [GIT].

Note in the plots below, time is defined in minutes.

The result is:

For the birth rate $\lambda$ and the death rate $\mu$, these will need to be approximated via simulations to see what mix of those parameters fit our expectations the best.

However for $a$, $C_{0}$, and $\gamma$, we can find values for those. According to the paper [Glioma], under the dose-response section of the Cytochalasin D graphics section, it was found that the IC50 value was $15.43$ micro-moles ($\mu$M units), meaning that the population would be cut in half after the 1 hour experiment. We can also observe this phenomena approximately from the graphic given in the paper [Glioma] included here:

Thus, we can conclude that $C_{0}=15.43\mu$M. We can adjust this parameter by changing $\beta=1$ (the treatment scaling factor) to $\beta=\frac{1}{100}$ to allow for the ODE model (2.2) to be interpretable in a graphical sense, so can be $C_{0}=0.1543\mu$M. Subsequently, it was also done for the simulations runs to keep our parameters consistent. Now for $a$, the initial population, the paper offers two initial populations, one of low density cell count of $1\times 10^{5}$ NPA cells, and one of high density cell count of $2\times 10^{5}$ NPA cells. At low density, no oncostreams were observed, whereas high density conditions facilitated robust oncostream formation. So it is safe to assume that our initial population $a=2\times 10^{5}$ NPA cells.

For $\gamma$, given the half-life time of Natural Cytochalasin D, as mentioned in the drug background information in section 2.2, we can calculate this decay rate. Given that the half life is 10 minutes we can calculate $\gamma$ as:

$$\displaystyle{\text{Given that: }N_{0}\left(\dfrac{1}{2}\right)^{\dfrac{t}{t_{\frac{1}{2}}}}=N_{0}e^{-\gamma t}\text{ are both exponential decay equations}}$$

$$\displaystyle{t_{\frac{1}{2}}\text{ represents half-life, we can solve for }\gamma}$$

(9)

$$\begin{split}\left(\dfrac{1}{2}\right)^{\dfrac{t}{t_{\frac{1}{2}}}}&=e^{-\gamma t}\\ \ln\left(\left(\dfrac{1}{2}\right)^{\dfrac{t}{t_{\frac{1}{2}}}}\right)&=\ln\left(e^{-\gamma t}\right)\\ \frac{t}{t_{\frac{1}{2}}}\ln\left(\frac{1}{2}\right)&=-\gamma t\\ -\frac{t}{t_{\frac{1}{2}}}\ln(2)&=-\gamma t\\ \frac{\ln(2)}{t_{\frac{1}{2}}}&=\gamma\\ \end{split}$$

From here, we can explicitly calculate $\gamma_{1}$ as:

$$\displaystyle{\gamma_{1}=\frac{\ln(2)}{10\text{ minutes}}\approx\frac{0.0693}{\text{minute}}}$$

The same process can be followed for Liposomal Cytochalasin Treatment to find that:

$$\displaystyle{\gamma_{2}=\frac{\ln(2)}{10\text{ minutes}}\approx\frac{0.0029}{\text{minute}}}$$

So to reiterate, our known parameters for Natural Cytochalasin Treatment are:

	$\displaystyle a$	$\displaystyle=2\times 10^{5}\text{ NPA cells}$
	$\displaystyle C_{0}$	$\displaystyle=0.1543\text{ micro-moles }(\mu\text{M})$
	$\displaystyle\gamma_{1}$	$\displaystyle=\frac{0.0693}{\text{minute}}$

And the for Liposomal Cytochalasin Treatment, the only difference is:

$\displaystyle\gamma_{2}$

$\displaystyle=\frac{0.0029}{\text{minute}}$

For this process, we can use a known relation and then solve for $f_{T}(t)$, the probability density function of the time, $T$, the time until the next event.

We know that:

$$\displaystyle{\mathbb{P}(T>t+\Delta t)=\mathbb{P}(T>t)\mathbb{P}_{\text{no event during $\Delta t$}}}$$

Which can be written as:

$$\displaystyle{\mathbb{P}(T>t+\Delta t)=\biggr{(}1-(\lambda+\mu+C(t))n\Delta t\biggl{)}\mathbb{P}(T>t)}$$

$$\displaystyle{\mathbb{P}(T>t+\Delta t)=\biggr{(}1-(\lambda+\mu+C_{0}e^{-\gamma t})n\Delta t\biggl{)}\mathbb{P}(T>t)}$$

$$\displaystyle{\mathbb{P}(T>t+\Delta t)-\mathbb{P}(T>t)=-(\lambda+\mu+C_{0}e^{-\gamma t})n\Delta t\mathbb{P}(T>t)}$$

$$\displaystyle{\frac{\mathbb{P}(T>t+\Delta t)-\mathbb{P}(T>t)}{\Delta t}=\frac{-(\lambda+\mu+C_{0}e^{-\gamma t})n\Delta t\mathbb{P}(T>t)}{\Delta t}}$$

This is the form of the derivative, thus:

$$\displaystyle{\frac{d{P}(T>t)}{dt}=-(\lambda+\mu+C_{0}e^{-\gamma t})n\mathbb{P}(T>t)}$$

And we can easily solve this ODE, with $A$ being an integration constant:

$$\displaystyle{\mathbb{P}(T>t)=A\exp\left[-(\lambda+\mu)nt+\dfrac{C_{0}ne^{-\gamma t}}{\gamma}\right],A\in\mathbb{R}}$$

We can solve for $A$ by using the fact that $\mathbb{P}(T>0)=1$:

$$\displaystyle{\mathbb{P}(T>0)=1=A\exp\left[0+\dfrac{C_{0}ne^{0}}{\gamma}\right]}$$

$$\displaystyle{1=A\exp\left[\dfrac{C_{0}n}{\gamma}\right]}$$

$$\displaystyle{e^{-\dfrac{C_{0}n}{\gamma}}=A}$$

Then we can express this as the cumulative distribution function of $T$:

$$\displaystyle{F_{T}(t)=1-\mathbb{P}(T>t)=1-e^{-\dfrac{C_{0}n}{\gamma}}\exp\left[-(\lambda+\mu)nt+\dfrac{C_{0}ne^{-\gamma t}}{\gamma}\right]}$$

$$\displaystyle{F_{T}(t)=1-\mathbb{P}(T>t)=1-\exp\left[-(\lambda+\mu)nt+\dfrac{C_{0}n(e^{-\gamma t}-1)}{\gamma}\right]}$$

And then we can get the probability density function by taking the derivative of $F_{T}(t)$:

$$\displaystyle{f_{T}(t)=n\left(C_{0}e^{-\gamma t}+\lambda+\mu\right)\exp\left[-(\lambda+\mu)nt+\frac{C_{0}n(e^{-\gamma t}-1)}{\gamma}\right]}$$

As a sanity check, we can confirm that the probability density function is legitimate by integrating it. Using Mathematica (code in Appendix 10) we can see that the integral produces 1, meaning it is a legitimate probability density function.

$$\displaystyle{\int_{0}^{\infty}f_{T}(t)dt=\int_{0}^{\infty}n\left(C_{0}e^{-\gamma t}+\lambda+\mu\right)\exp\left[-(\lambda+\mu)nt+\frac{C_{0}n(e^{-\gamma t}-1)}{\gamma}\right]dt=1}$$

Here are graphs representing the distribution for both Natural and Liposomal Cytochalasin D. The y-axis was scaled down by a factor of 100 to show probabilities a decimal to better represent what the probability is. A population of $10$ was used to allow better representation of and visualization of the PDF curve. Both of these plots were generated using code that can be found in Appendix 12.

From these plots, it is easy to observe that the probability of the next event happening is highest when $0<\Delta t<<1$. As a result of this, births and deaths will occur extremely frequently with very little time in between each birth or death.

Simulation results of 10 different runs yielded these graphs (using MATLAB code found in Appendix 12).

In this Birth-Death simulation with Liposomal Cytochalasin D Treatment, we can see that on average every simulation run has the same downward trend and actually reduces the population size as time passes. This is very good because this indicates that the treatment is working properly and shows promising results for controlling and reducing the Cell Population.

In this Birth-Death Process simulation with Natural Cytochalasin D Treatment, we can see that it doesn’t actually work as much as the previous type of treatment. With the Population growing exponentially in some runs, and the overall trend still showing that the population grows.