September 25, 2017, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Modeling Stem/ Progenitor Cell-induced Neovascularization and Oxygenation Around Solid Implants

From MathBio

Jump to: navigation, search

Student recreation/ extension and article review by Matthew Cullen, Caitlin Lowe, Sam Verplanck,

Article "Modeling Stem/ Progenitor Cell-Induced Neovascularization And Oxygenation Around Solid Implants" by Harsh Vardhan Jain, Nicanor I Modovan, Helen M. Byrne


The author's goal was to create a mathematical model to explain the vascularization and oxygen tension found in a hydrogel device 1. The authors did not do the experiment themselves, but created a model that fit the data from Butt et al 2. The authors use a series of reaction-diffusion-aversion reactions to create their model. This model explores how the device treated with stem cells creates an adipogenesis reaction and without stem cells a foreign body response (FBR) is found. Adipogenesis is the process of pre-adipocytes becoming adipocytes through cell differentiation. This adipogenesis is assumed to indicate neovascularization, the creation of blood vessels. FBR is a wound healing response where fibroblasts and inflammatory cells are formed. The FBR response caused a lower oxygen content to the device and was verified by the model. The adipogenisis response allowed a higher oxygen tension and the tension increased with the thickness of stem cells. The change in tension was negligible with a higher density of stem cells. The results did show that after a certain thickness the vascularization takes too long and is considered a failure.

Context / Biological phenomenon under consideration

Medical devices such as tissue engineering constructs and solid implants are often rejected by the host. This rejection is caused by the foreign body reaction (FBR). There is also an avascular sheath formed, which limits the amount of molecular oxygen to reach the device. The biomedical devices need oxygen to operate properly. In this paper a model was created based on the paper "Stimulation of peri-implant vascularization with bone marrow-derived progenitor cells: monitoring by in vivo EPR oximetry" 2. The model was created to give a possible mathematical explaination to the research. In this model, two sub models are discussed. The FBR submodel deals with the inflammatory cell density. The adipogenesis model analyses how the density of stem cells and pre-adipocyte density affect the vascularity and oxygen tension in the device.


Stem cells and progenitor cells have been a popular and controversial topic over the past view decades. There have been a large number of scientific papers regarding stem cell research and applications. One of the most important preceding papers to this paper is regarding simulation of peri-implant vasularization by Butt et al. 2. In this paper the scientists ran an experiment by setting up a oxygen sensing probe inside an implant. Then the researchers surrounded that implant with a layer of marrow-derived progenitor cells. After the implant, the researchers monitored oxygen pressure over a time span of ten weeks. A main goal of Jain's paper is to develop a mathematical model that mimics the experimental results from Butt's paper. This paper allowed Jain et al. to estimate values for key parameters in their model. Also, all simulations were ran over a 10-week timescale corresponding to the timescale of the experiments in Butt et al.

The mathematical model discussed in paper was based on several models of wound healing angiogenesis, 3 4 5 6 and glucose transport to implanted glucose sensors, 7 8. The earliest of these papers date back to 1996 while the most recent were published in 2010.

Mathematical Model

The model is of peri-implant vascularization, and provides a possible explanation for experimental observations. This model is comprised of ordinary and partial differential equations, that create a system of Reaction diffusion equations. The model is divided into three sections Core Model Equations, FBR submodel, and adipogenesis submodel.

  1. Core Model Equations

In the model the implant is cylindical with oxygen sensitive crystals located at one end. The geometry is considered to be two-dimensional and radially symmetric with a radius of "r".

The equation that governs the microvascular density V(r,t) comes from combining the random movement of endothelial cells, chemotaxis, formation, and remodeling of vessel cells:

{\frac  {\delta V}{\delta t}}={\frac  {D_{v}}{r}}{\frac  {\delta }{\delta r}}\left(r{\frac  {\delta V}{\delta r}}\right)-{\frac  {1}{r}}{\frac  {\delta }{\delta r}}\left(r\;\chi (U)V{\frac  {\delta U}{\delta r}}\right)+f_{0}(U,V)-\rho V^{2} (1)


The chemotactic sensitivity is assumed to decrease as angiogenic factors increase and is given by given by receptor kinetic laws: \chi (U)={\frac  {\chi _{o}}{(U+K_{u})^{2}}}

The vascular formation rate is assumed to be increasing and the function becomes:

f_{v}(U,V)=\rho \left[V_{h}+{\frac  {\sigma _{1}U}{\sigma _{2}+U}}\right]V

The vascular remodeling term is in equation (1)\rho V^{2} is needed so that the steady state U=0relates to the metabolic demands of healthy tissue V=V_{h}

Assuming that the oxygen inside the tissue is delivered by vascular network and diffusion, and then consumed by the tissue the oxygen tension C(r,t) becomes:

{\frac  {\delta C}{\delta t}}={\frac  {D_{c}}{r}}{\frac  {\delta }{\delta r}}\left(r{\frac  {\delta C}{\delta r}}\right)+f_{c}(V,C)-\lambda _{c}C (2)


The rate that the oxygen is delivered by the vascular system is: f_{c}(V,C)=\mu V(C_{v}-C)r This rate is assumed to be proportional to the vessel density V, and the difference between oxygen tension in the vasculature and tissue.

The Angiogenic factor U(r,t) is given by :

{\frac  {\delta U}{\delta t}}={\frac  {D_{u}}{r}}{\frac  {\delta }{\delta r}}\left(r{\frac  {\delta U}{\delta r}}\right)+f_{u}(.)-\lambda _{u}U-\beta VU (3)


For this equation it is assumed that the angiogenic factors are diffused and produced by the tissue, and then removed by the vasculature.

Initial and Boundary Conditions

Because the hermetically oxygen sensor is sealed the flux at the implant edge is zero: {\frac  {\delta V}{\delta r}}={\frac  {\delta C}{\delta r}}={\frac  {\delta U}{\delta r}}=0 when evaluated at r=R_{i}

The assumption was made that V and C approach normal values and U goes to 0: V(\infty )=V_{h},\;C(\infty )=C_{h},\;U(\infty )=0

Finally, it is assumed that at t=0 inside of the device there are no vessels or oxygen, and the tissue does not have any angiogenic factors.


  1. FBR submodel

The vascularization is assumed to happen in two parts. The first phase is the invasion of the hydrogel by inflammatory cells. M(r,t) is the inflammatory cell density which is given by the invasion in to the hydrogel plug equal to the recruitment and differentiation in implant site minus the natural death of inflammatory cells.

{\frac  {\delta M}{\delta t}}+H(t_{f}-t)\lambda {\frac  {\delta M}{\delta r}}=H(t-t_{f})\left({\frac  {a_{m}}{(b_{m}+G)}}M-\lambda _{m}M\right) (4)


The second part of the vascularization comes from the avascular sheath formation. G(t) is the avascular sheath density. The formation of this sheath happens because the implant is made of hard plastic, stimulates the FBR reaction. The formation is given by: {\frac  {dG}{dt}}=H(t-t_{f})c_{m}M (5)

The H(.) function found in both equation (4) and (5) is the Heaviside step function: H(t-a)=1,\;t>a and H(t-a)=0,\;t<a

The control experiments were done without stem cell supplementation. The FBR submodel is joined to the core model by the angiogenic facter production term : f_{f}(M,t,r)=K_{m}M (6)

The boundary and initial conditions are:


M(r,0)=0,\;r\neq L_{0}

M(\infty ,t)=0



  1. Adipogenesis submodel

In the Adipogenesis submodel, the stem cells are assumed to be uniformly distributed in the hydrogel plug. The cells are also assumed to stimulate an adipogenic response rather than FBR. The stem cell density S(r,t) and pre-adipocyte density A(r,t) are defined on R_{i}<r<\infty and governed by equations (7) and (8).

Equation (7) represents the death of the stem cells:

{\frac  {\delta S}{\delta t}}=-\lambda _{s}S(7)

Equation (8) is the Invasion into hydrogel plug equal to the Proliferation and Differentiation of the stem cells:

{\frac  {\delta A}{\delta t}}+H(t_{a}-t)\left({\frac  {\lambda _{A}S}{K_{A}+S}}\right){\frac  {\delta }{\delta r}}(A)={\frac  {a_{A}S}{b_{a}+S}}A-c_{A}A(8)

There is a lack of experimental data so the assumption is made that equation (8) is a direct response to the presence of stem cells.

The assumption was also made that the pre-adipocytes proliferate because of stem cells at an increasing rate.

The adipogenesis submodel joins the core model in equation 3 by replacing f_{u} with: f_{A}(A,t,r)=\xi _{A}A

The initial and boundary conditions are given by:




A(r,0)=0,\;r\neq L_{0}


Example Discussion

Oxygen tension C(r,t) is one of the most important factors in this model. Let's discuss what the PDE describing oxygen tension actually means.

What does the {\frac  {D_{c}}{r}}*{\frac  {\delta }{\delta r}}\left(r*{\frac  {\delta C}{\delta r}}\right) term represent?

This term represents diffusion of oxygen. Note that the first time derivative is set equal to the second spatial derivative which is very common PDE that also shows how heat is transferred in 1 dimensional cylindrical coordinates. The term D_{c} is the diffusion coefficient; a higher diffusion coefficient is indicative of a larger rate of diffusion. The oxygen diffusion coefficient used in this model is 207.36 mm2 per day.

What does the f_{c}(V,C)=\mu V(C_{v}-C)r function represent?

This function describes the delivery rate of oxygen from the vasculature to the tissue. In this equation \mu is 14 per microvessel density per day. A larger \mu corresponds to a higher delivery rate of oxygen. Here, Cv represents the microvascular oxygen tension. The value used is 20 mmHg. Note, this term is in summation with diffusion because a higher delivery rate corresponds with a higher rate of change in oxygen tension.

What does the -\lambda _{c}C term represent?

This term represents the consumption of oxygen by the tissue. Here, \lambda _{c} represents the tissue uptake rate of oxygen. The value used was 42 per day. Note, this term is negative because a higher uptake rate would lead to a smaller rate of change of oxygen tension.

Published Results

In order, to make predictions Jain et al. first had to fit the parameters the models to experimental data from Butt et al. 2. The experimental data that was fitted is shown in Figure 1 below.

Oxy press.Jpg

Fig. 1: Oxygen Pressure vs. Time of (A) the control experiments and (B) the cell-treatment experiments

Once the appropriate values for key parameters were chosen, Jain et al. produced predictions of the effect of varying the thickness of stem cell coating (scaffold).

Scaffold thick.Jpg

Fig. 2: Predicted behavior of (A) maximum oxygen tension at the implant edge (B) time elapsed to reach maximum oxygen tension and (C) the oxygen tension at the implant edge after 10 weeks of implantation vs. scaffold thickness

Jain et al. ran more simulations with scaffold thickness held constant at 1.6 mm, while varying stem cell density. The results are shown below.

Cell density.Jpg

Fig. 3: Predicted behavior of (A) maximum oxygen tension at the implant edge (B) time elapsed to reach maximum oxygen tension and (C) the oxygen tension at the implant edge after 10 weeks of implantation vs. stem cell density

The final important result from this paper is in understanding the limitations of a stem cell approach. A key value in this study is ta, the time it takes for the leading edge of the pre-adipocytes to reach the implant. If this value approaches infinity than any treatment would be seen as a failure. Shown below are how ta vary with thickness of stem cell coating and density of stem cells.


Fig. 4: Time for invading pre-adipocytes to reach implanted device edge vs. (A) scaffold thickness and (B) stem cell density

Student Produced Results

A team of three students enrolled in APPM 4390 at the University of Colorado attempted to reproduce results from the paper. The main function used to find these results was pdepe in MATLAB. See the appendix for code. Below are a few figures generated from the FBR submodel. FBR M.jpg

Fig. 5: Inflammatory Cell Density as a function of time (days) and distance (mm)

FBR V.jpg

Fig. 6: Vascularization as a function of time (days) and distance (mm)

Note in Figure 5 as r approaches infinity the inflammatory cell density approaches 0. This is consistent with the boundary condition and physical system. Figure 6 shows vascularization as a function of space and time. Note, vascularization only occurs in the Mitrogel layer surrounding the implant. That is, from radii ranging from 2.0 mm to 3.6 mm

FBR C.jpg

Fig. 7: Oxygen Tension as a function of time (days) and distance (mm)

FBR U.jpg

Fig. 8: Angiogenic Factor Concentration as a function of time (days) and distance (mm)

Shown above in Figure 7 shows how oxygen tension depends on time and distance. Note as r increases, oxygen tension approaches the body's healthy level. Also, as time increases, the oxygen tension increases asymptotically. Figure 8 shows angiogenic factor concentration. As the distance from the implant increases, U approaches zero. This result makes sense for a healthy organism. Also, there is a wave in angiogenic factor concentration with its peak around day 35.

FBR 3t.jpg

Fig. 9: Vascularization, Oxygen Tension, and Angiogenic Factor Concentration as a function of time (days) at the edge of the oxygen sensor

FBR 3r.jpg

Fig. 10: Vascularization, Oxygen Tension, and Angiogenic Factor Concentration as a function of distance from center of probe (mm) at 5 weeks.

Figures 9 and 10 show 2D representations of the 3D meshplots above. Note that the middle graph in Figure 9 is analogous to Figure 1A. Both representations show an increase vascularization to begin. However, the magnitudes do not correspond well. This could be due to the student's error in inputting values, although the code was thoroughly checked.

After producing a representation of the FBR submodel, the team of students attempted to explore the Adipogenesis submodel. After inputting the the set of partial differential equations described in the Mathematical Model section of the paper into MATLAB and defining parameters, the output seemed to be numerically unstable. Below, Figure 11 shows values outputted for the pre-adipocyte density.


Fig. 11:A misrepresentation of pre-adipocyte density as a function of time (days) and distance (mm)

In addition to an unrealistic shape, the plot above has negative density values. Thus it is a clearly not a valid representation.

Analysis/Discussion of Published Results

The effect of thickness of the stem cell coating is quite interesting. A thicker coating implies more stem cell so there is an increase in maximum oxygen tension shown in Figure 2A. However, a maximum is reached around 2.2 mm. Past this point oxygen tension falls with an increased thickness. Figure 2B shows the 2.2 mm thickness threshold in a different light. Here the time taken to reach the maximum oxygen tension decreases with an increasing thickness until that 2.2 mm thickness threshold. Figure 3C shows an increasing oxygen tension after ten weeks with an increasing thickness.

Stem cell density is also an important variable to look at. Shown in Figure 3A is how increasing stem cell density does not greatly increase the maximum oxygen tension. On the other, increasing stem cell density can greatly increase the time elapsed until that maximum value is reached (Fig. 3B). The increase in time taken to achieve a maximum oxygen tension corresponds to Figure 4C where there is a higher oxygen tension after 10 weeks with an increasing density.

A very useful result from this model is show in Figure 4. That is, there are distinct limitations to the use of stem cells. The model reaches an asymptote at a critical thickness of stem cells. At this thickness the pre-adipocytes take too long to reach the implant and that is seen as a failure. A similar phenomenon occurs at a minimum stem cell density. Below this density there are not enough stem cells in an area to allow for pre-adipocytes to travel to the implant.

More recently Butt et al have developed a new oxygen tension measuring sensor. The mathematical model was applied to the new results and "excellent qualitative and and quantitative results has been obtained"1. This agreement implies that the biological mechanics are described well with the mathematical model.

Analysis/Discussion of Student Results

The FBR submodel's reproduced results follow the published results fairly well. All boundary conditions are satisfied and general behavior is similar to the published results. The main difference is in magnitude of the results. Many numerical values differ from those published in Jain et al. One cause could be input error of the many parameters. The team of students checked these values many times and they appear to be correct. Another cause could be a mistyped value in the scientific report.

The team was unable to reproduce valid results for the adipogenesis submodel. A dead give away for this was in negative pre-adipocyte density values. One cause of this could be the set of PDE's is numerically unstable when evaluating solutions using MATLAB's pdepe function.

Conclusion and Further Study

The team of students was able tor reproduce consistent parts of the FBR submodel defined by Jain et al. The adipogenesis submodel showed to be numerically unstable with solver used. Everyone on the World Wide Web is invited and encouraged to pick through the code used. It can be found in the appendix. Further study could be done by optimizing code, using different numerical solvers, and performing a sensitivity analysis.

External Links

1. Jain, Harsh Vardhan, Nicanor I. Moldovan, and Helen M. Byrne. "Modeling Stem/Progenitor Cell-Induced Neovascularization and Oxygenation Around Solid Implants." Tissue Engineering Part C: Methods 18.7 (2012): 487-495. 1

2. Butt, O.I., Carruth, R., Kutala, V.K., Kuppusamy, P., and Moldovan, N.I. Stimulation of peri-implant vascularization with bone marrow-derived progenitor cells: monitoring by in vivo EPR oximetry. Tissue Eng 13, 2053, 2007. 2

3. Pettet, G.J., Byrne, H.M., McElwain, D.L., and Norbury, J. A model of wound-healing angiogenesis in soft tissue. Math Biosci 136, 35, 1996. 3

4. Pettet, G.J., Chaplain, M.A., McElwain, D.L., and Byrne, H.M. On the role of angiogenesis in wound healing. Proc Biol Sci 22, 263, 1996. 4

5. Schugart, R.C., Friedman, A., Zhao, R., and Sen, C.K. Wound angiogenesis as a function of tissue oxygen tension: a mathematical model. Proc Natl Acad Sci U S A 105, 2628, 2008. 5

6. Stamper, I.J., Byrne, H.M., Owen, M.R., and Maini, P.K. Modelling the role of angiogenesis and vasculogenesis in solid tumour growth. Bull Math Biol 69, 2737, 2007. 6

7. Novak, M.T., Yuan, F., and Reichert, W.M. Modeling the relative impact of capsular tissue effects on implanted glucose sensor time lag and signal attenuation. Anal Bioanal Chem 398, 1695, 2010. 7

8. Singh, S., and McShane, M. Enhancing the longevity of microparticle-based glucose sensors towards one month continuous operation. Biosens Bioelectron 25, 1075, 2010. 8


The following code was used to produce the student results

function VCU_calc
%This function will solve the system of PDEs with the parameters from the 
%original paper

close all

%Set up domain and important constants
m = 1;
r = 2:.05:25; % from Ri to Lo in mm
t = 0:.1:80;% days
global M dt dr
dt = (t(end)-t(1))/length(t);
dr = (r(end)-r(1))/length(r);
M = Mfind(t,r);

%Solve the PDEs
sol = pdepe(m,@pdex9pde,@pdex9ic,@pdex9bc,r,t);

%Extract solutions
V = sol(:,:,1);
C = sol(:,:,2);
U = sol(:,:,3);

%Plot relevant information
title('V (r,t)')
xlabel('Distance r')
ylabel('Time t')

title('C (r,t)')
xlabel('Distance r')
ylabel('Time t')

title('U (r,t)')
xlabel('Distance r')
ylabel('Time t')

xlabel('Time t')

xlabel('Time t')

xlabel('Time t')

title('V(r,5 weeks)')
xlabel('distance r')

title('C(r,5 weeks)')
xlabel('distance r')

title('U(r,5 weeks)')
xlabel('distance r')

% --------------------------------------------------------------
function [c,f,s] = pdex9pde(r,t,u,DuDr)

global M dt dr
D_v = .009; % mm^2/day
chi_0 = 2.5; % mm^2/day
K_U = 2; % micrometers
D_c = 207.36; % mm^2/day
D_u = 6.048; % mm^2/day
rho = .1; % per microvessel density per day
V_h = 1.035; % number (in thousands) mm^2
sigma1 = 57; % microvessel density per day
sigma2 = 2; % micrometers
mu = 14; % per microvessel density per day
c_v = 20; % mmHg
lambda_c = 42; % per day 
k_m = 1; % microM per inflammatory cell density per day
lambda_u = 15.547; % per day
beta = .1; % per microvessel density per day
chi = chi_0/((u(3)+K_U).^2);

%Find the needed index of t and r to find the closest M value
if t/dt < 801
    tindex = floor(t/dt + 1);
    tindex = 801;

f_fbr = k_m*M(tindex,round((r-2)/dr)+1);

%The PDE equations

c = [1; 1;1]; 

f = [DuDr(1).*D_v - chi.*u(1)*DuDr(3);D_c.*DuDr(2);D_u.*DuDr(3)];

s = [rho*(V_h + sigma1*u(3)/(sigma2+u(3))*u(1) - u(1)^2) ...
    mu.*u(1).*(c_v-u(2)).*r - lambda_c.*u(2)...

% --------------------------------------------------------------

function u0 = pdex9ic(r)
%Initial conditions

V_h = 1.035; % number (in thousands) mm^2
C_h = 5.130; % mmHg

u0 = [heaviside(r-3.6).*V_h; heaviside(r-3.6)*C_h; 0]; 

% --------------------------------------------------------------

function [pl,ql,pr,qr] = pdex9bc(~,~,~,ur,~)
%Boundary conditions

D_v = .009; % mm^2/day
D_c = 207.36; % mm^2/day
D_u = 6.048; % mm^2/day
V_h = 1.035; % number (in thousands) mm^2
C_h = 5.130; % mmHg

pl = [0;0;0]; 
ql = [1/D_v;1/D_c; 1/D_u]; 
pr = [ur(1)-V_h; ur(2) - C_h; ur(3)]; 
qr = [0; 0; 0]; 


function M = Mfind(t,r)
%This function will solve the PDE for M over a given domain.

lamdaF = .237;
M0 = .07;
L0 = 3.6;
tf = 7;
dt = (t(end)-t(1))/length(t);
dr = (r(end)-r(1))/length(r);
M = zeros(length(t),length(r));

%solve until tf
tfindex = floor(tf/dt);

for ii = 1 : tfindex    
    index = ceil((L0 - t(ii)*lamdaF)/dr);
    M(ii,index-39) = M0;

M(tfindex+ 1,:) = M(tfindex,:);

%Solve after tf
[~, MGstuff] = ode45(@MvsG, t(tfindex+2:end),[M(tfindex+1,1) 0]');

M(tfindex+2:end,1) = MGstuff(:,1);

%Plot the results
xlabel('Distance r')
ylabel('Time t')

function du = MvsG(~,u)
%coupled ODEs of M and G

du = [0 0]';

aM = .024;
bM = .061;
lamdaM = .237;
cM = .007;

du(1) = ((aM/(bM + u(2))) - lamdaM) *u(1);
du(2) = cM * u(1);


function Afind
%This function solves for the pre-adipocyte density. It was found to yield 
%non realistic answers due to numeric instability.

close all

r = 2:.05:50; % from Ri to Lo in mm
t = 0:.1:80;% days

L0 = 3.6;
S0 = 10;
lamdaS = .099;
m = 0;
ta = 7;
aA = .414;
bA = 3.792;
cA = .063;

global S dt dr

dt = (t(end)-t(1))/length(t);
dr = (r(end)-r(1))/length(r);
taindex = floor(ta/dt);
S = zeros(length(t),length(r));
A = S;

%assign all values of S
for ii = 1:length(t)    
    S(ii,1:(floor(L0/dr))) = S0*exp(-lamdaS*t(ii));

%Solve for A until ta 
A(1:taindex,:) = pdepe(m,@pdex9pde,@pdex9ic,@pdex9bc,r,t(1:taindex));

%Solve for the remaining values of A
for ii = taindex+1:length(t)
    for jj = 1:length(r)
        S1 = S(ii, jj);
        A(ii,jj) = A(taindex,jj)*exp(t(ii-taindex)*((S1*aA/(bA + S1) - cA)));

xlabel('Distance r')
ylabel('Time t')


function [c,f,s] = pdex9pde(r,t,u,~)
global S dt dr
Ka = 5;
lamdaa = .130;
aA = .414;
bA = 3.792;
cA = .063;

%Determine needed index of t
ta = 7;
taindex = floor(ta/dt);

if t < taindex
    tindex = floor(t/dt)+1;
    tindex = taindex;

%Determine needed value from S
S1 = S(tindex, round(r/dr-40));

c = 1;
f = -((lamdaa*S1)/(Ka + S1))*u;
s = (S1*aA/(bA + S1) - cA)*u(1);


function u0 = pdex9ic(r)
%Initial conditions

global dr
L0 = 3.6;
A0 = .35;

if floor((r-L0)/dr) == 0
    u0 = A0;
    u0 = 0;


function [pl,ql,pr,qr] = pdex9bc(~,ul,~,ur,~)
%boundary conditions
pl = ul;
ql = 0;
pr = ur;
qr = 0;