
MBW:Modeling Stem/ Progenitor Cellinduced Neovascularization and Oxygenation Around Solid ImplantsFrom MathBioStudent recreation/ extension and article review by Matthew Cullen, Caitlin Lowe, Sam Verplanck, Article "Modeling Stem/ Progenitor CellInduced Neovascularization And Oxygenation Around Solid Implants" by Harsh Vardhan Jain, Nicanor I Modovan, Helen M. Byrne Contents
SummaryThe 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 reactiondiffusionaversion 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 preadipocytes 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 considerationMedical 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 periimplant vascularization with bone marrowderived 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 preadipocyte density affect the vascularity and oxygen tension in the device. HistoryStem 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 periimplant 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 marrowderived 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 10week 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 ModelThe model is of periimplant 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.
In the model the implant is cylindical with oxygen sensitive crystals located at one end. The geometry is considered to be twodimensional 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: (1)
The chemotactic sensitivity is assumed to decrease as angiogenic factors increase and is given by given by receptor kinetic laws: The vascular formation rate is assumed to be increasing and the function becomes:
The vascular remodeling term is in equation (1) is needed so that the steady state relates to the metabolic demands of healthy tissue 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: (2)
The rate that the oxygen is delivered by the vascular system is: 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 : (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: when evaluated at The assumption was made that V and C approach normal values and U goes to 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.
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. (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: (5) The H(.) function found in both equation (4) and (5) is the Heaviside step function: and The control experiments were done without stem cell supplementation. The FBR submodel is joined to the core model by the angiogenic facter production term : (6) The boundary and initial conditions are:
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 preadipocyte density A(r,t) are defined on and governed by equations (7) and (8). Equation (7) represents the death of the stem cells: (7) Equation (8) is the Invasion into hydrogel plug equal to the Proliferation and Differentiation of the stem cells: (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 preadipocytes proliferate because of stem cells at an increasing rate. The adipogenesis submodel joins the core model in equation 3 by replacing with: The initial and boundary conditions are given by:
Example DiscussionOxygen 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 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 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 mm^{2} per day.
What does the function represent? This function describes the delivery rate of oxygen from the vasculature to the tissue. In this equation is 14 per microvessel density per day. A larger corresponds to a higher delivery rate of oxygen. Here, C_{v} 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 term represent? This term represents the consumption of oxygen by the tissue. Here, 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 ResultsIn 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. Fig. 1: Oxygen Pressure vs. Time of (A) the control experiments and (B) the celltreatment 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). 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. 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 t_{a}, the time it takes for the leading edge of the preadipocytes to reach the implant. If this value approaches infinity than any treatment would be seen as a failure. Shown below are how t_{a} vary with thickness of stem cell coating and density of stem cells. Fig. 4: Time for invading preadipocytes to reach implanted device edge vs. (A) scaffold thickness and (B) stem cell density Student Produced ResultsA 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. Fig. 5: Inflammatory Cell Density as a function of time (days) and distance (mm) 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 Fig. 7: Oxygen Tension as a function of time (days) and distance (mm) 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. Fig. 9: Vascularization, Oxygen Tension, and Angiogenic Factor Concentration as a function of time (days) at the edge of the oxygen sensor 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 preadipocyte density. Fig. 11:A misrepresentation of preadipocyte 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 ResultsThe 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 preadipocytes 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 preadipocytes 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 ResultsThe 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 preadipocyte 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 StudyThe 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 Links1. Jain, Harsh Vardhan, Nicanor I. Moldovan, and Helen M. Byrne. "Modeling Stem/Progenitor CellInduced Neovascularization and Oxygenation Around Solid Implants." Tissue Engineering Part C: Methods 18.7 (2012): 487495. 1 2. Butt, O.I., Carruth, R., Kutala, V.K., Kuppusamy, P., and Moldovan, N.I. Stimulation of periimplant vascularization with bone marrowderived 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 woundhealing 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 microparticlebased glucose sensors towards one month continuous operation. Biosens Bioelectron 25, 1075, 2010. 8 AppendixThe 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 figure mesh(r,t,V) title('V (r,t)') xlabel('Distance r') ylabel('Time t') figure mesh(r,t,C) title('C (r,t)') xlabel('Distance r') ylabel('Time t') figure mesh(r,t,U) title('U (r,t)') xlabel('Distance r') ylabel('Time t') figure subplot(3,1,1) plot(t,V(:,1)) title('V(Ri,t)') xlabel('Time t') subplot(3,1,2) plot(t,C(:,1)) title('C(Ri,t)') xlabel('Time t') subplot(3,1,3) plot(t,U(:,1)) title('U(Ri,t)') xlabel('Time t') figure subplot(3,1,1) plot(r,V(700,:)) title('V(r,5 weeks)') xlabel('distance r') subplot(3,1,2) plot(r,C(700,:)) title('C(r,5 weeks)') xlabel('distance r') subplot(3,1,3) plot(r,U(350,:)) title('U(r,5 weeks)') xlabel('distance r') %  function [c,f,s] = pdex9pde(r,t,u,DuDr) %Constants 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); else tindex = 801; end f_fbr = k_m*M(tindex,round((r2)/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_vu(2)).*r  lambda_c.*u(2)... f_fbrlambda_u*u(3)beta.*u(1).*u(3)]'; %  function u0 = pdex9ic(r) %Initial conditions V_h = 1.035; % number (in thousands) mm^2 C_h = 5.130; % mmHg u0 = [heaviside(r3.6).*V_h; heaviside(r3.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. %Constants 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,index39) = M0; end 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 mesh(r,t,M) title('M(r,t)') 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 preadipocyte 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 %Constants 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)); end %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(iitaindex)*((S1*aA/(bA + S1)  cA))); end end mesh(r(1:60),t,A(:,1:60)) title('A(r,t)') xlabel('Distance r') ylabel('Time t') % function [c,f,s] = pdex9pde(r,t,u,~) %Constants 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; else tindex = taindex; end %Determine needed value from S S1 = S(tindex, round(r/dr40)); 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((rL0)/dr) == 0 u0 = A0; else u0 = 0; end % function [pl,ql,pr,qr] = pdex9bc(~,ul,~,ur,~) %boundary conditions pl = ul; ql = 0; pr = ur; qr = 0; 