In this paper, we present a dynamic Bayesian network (DBN) approach to modeling vascularization in engineered tissues. Injuries and diseases can cause significant tissue loss to the degree where the body is unable to heal itself. Tissue engineering aims to replace the lost tissue through use of stem cells and biomaterials. For tissue cells to multiply and migrate, they need to be close to blood vessels, and hence proper vascularization of the tissue is an essential component of the engineering process. We model vascularization through a DBN whose structure and parameters are elicited from experts. The DBN provides spatial and temporal probabilistic reasoning, enabling tissue engineers to test sensitivity of vascularization to various factors and gain useful insights into the vascularization process. We present initial results in this paper and then discuss a number of future research problems and challenges.
INTRODUCTION
People lose tissue due to accidents, medical
operations, treatments, and illnesses. While some organs,
e.g. liver, can replace the lost tissue most cannot
especially when the damage is too severe. For these kinds
of tissue damages, the lost tissue can be replaced by
engineering a new tissue through stem cells and
biomaterials [
An essential process for engineering a healthy tissue is the proper vascularization (formation of new blood vessels) of the tissue, as the tissue cells need to be close to the blood vessels both to discharge their waste and to receive nutrition and oxygen. The blood vessels need to spread out in the tissue, invade into the
The formation of new blood vessels are triggered and a↵ ected by growth factors that are released by distressed cells that are far from the existing blood vessels. When these growth factors reach existing blood vessels, they sprout new branches and these branches “search” for the distressed cells by following the gradient of the growth factor. This process, however, is stochastic for at least two reasons: i) even though growth factors are the main ingredients for causing sprouts, they are not the only elements that a↵ ect vascularization, and ii) the growth factors are increasingly more uniformly distributed as they go further away from the distressed cells, and hence the gradient is almost uniform, hindering the capability of the blood vessel finding its way correctly.
This inherent stochasticity in the vascularization process, the spatial nature of the tissue, and the temporal aspect of the vascularization make temporal graphical models a great fit for reasoning with uncertainty in vascularization. In this paper, we present a dynamic Bayesian network (DBN) for modeling vascularization in engineered tissues. We elicit the structure of the DBN from tissue engineering experts and we experiment with various parameter settings to provide further insights into the vascularization process. Because this is a first and novel application of DBNs to tissue engineering, it avails itself to many interesting future research directions and challenges.
Our contributions in this paper include: • We present a novel application of DBNs to vascularization in engineered tissues • We present initial results and insights, where we experiment with various parameter settings, and • We discuss several future research challenges and opportunities in detail.
The rest of the paper is organized as follows: in Section 2, we provide a brief background on tissue engineering and vascularization. In Section 3, we describe our DBN model for vascularization. We present our experimental setup and results in Section 4. In Section 5, we briefly discuss related work. We then discuss future research directions and challenges in detail in Section 6, and then conclude. 2
BACKGROUND In this section, we first provide a brief background on tissue engineering and vascularization and then discuss briefly why dynamic Bayesian networks (DBNs) are a good fit for modeling vascularization.
People lose tissue due to accidents, treatments, and illnesses. Some organs, e.g. liver, can replace the lost tissue while others cannot. Sometimes, the damage can be so severe that the body cannot heal itself. For example, bones can heal after smooth fractures. Yet, some fractures damage bone body so severely that the bone cannot regenerate. For these kinds of damages, the lost tissue can be replaced by engineering a new tissue through stem cells and biomaterials.
Stem cells are generic types of cells that have the ability to replicate and transform to any tissue. Stem cells, like all other cells, need to be close enough to the blood vessels so that they can forward their biological wastes to the vessels and they can be fed with nutrition and oxygen carried by the blood vessels. When a tissue is engineered through replication and transformation of stem and tissue cells, there is no existing blood vessel web in the environment; the only blood vessels available are the original vessels located at the edges, ready to sprout and progress to the depths of the newly-formed tissue.
The stem cells that do not have access to blood vessels will not be able to discharge waste and receive nutrition and oxygen. In such cases, a cell starts signaling about its needs by means of emitting chemicals called vascular endothelial growth factor (VEGF). VEGF diffuses and disperses in the environment. When it contacts a blood vessel, it triggers a new sprout of blood vessel towards the source of emission. The tip of these new sprouts typically follow the gradient of the VEGF to find the distressed cell. During this process, the newly-formed blood vessel can also branch and sprout new blood vessels. When the branches meet with other branches, they merge (this process is called anastomosis) and a blood circulation through the new vessel starts. The blood circulation helps nearby stem and tissue cells, which then stop emitting growth factors. This event is called angiogenesis or vascularization. Please see Figure 1 for an illustration of this process. Vascularization is a key process in tissue development. When cells that are emitting VEGF cannot be reached in time by the new blood vessels, the cells first fall in hypoxia (i.e., lack of Oxygen) and then start dying. Hence the formation of healthy tissue depends on appropriate vascularization; the blood vessels need to spread out in the newly-formed tissue, invade into the depth, and need to form connections to allow blood circulation.
Though it is well-known that the VEGF is a major contributor to sprouting of new blood vessels and that the tip of the blood vessel typically follows the gradient of the VEGF, there are still unknown factors that affect vascularization. Moreover, the VEGF distribution becomes more uniform as we get further away from the source of the emission and hence the gradient does not necessarily point to the distressed cell. Therefore, given our knowledge of the VEGF distribution the environment, the blood vessels do not necessarily follow a deterministic path; they also do a bit of exploration. This is where the uncertainty reasoning capabilities of probabilistic graphical models become handy for modeling vascularization.
In this paper, we model the vascularization process through dynamic Bayesian networks (DBNs) to enable tissue engineering researchers to reason with spatial and temporal growth of blood vessels. With the help of DBNs, the researchers can formulate and query the DBNs and try a number of parameter settings, without the need to experiment with every one of them in the lab. This process allows the researchers to gain further insights and formulate new in-vivo (on animals) and in-vitro (on glass) experiments. In this section, we describe our DBN model for vascularization.
We made a number of assumptions to
simplify the model. In this model, we assume a 2D
structure, whereas in real-life scenarios, the tissue
obviously has a 3D structure. In this 2D structure, which
is illustrated in Figure 2, as also assumed in [
Therefore, the status of a location at time t depends on: i) its status at time t, and ii) the statuses of its south neighbors at t. t ( ) t+1 ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( ) location Ltxy, we will drop the subscripts and hence simply use Lt, and when we refer to its neighbors t t t at its south L(x 1)(y 1), Lx(y 1), and L(x+1)(y 1) we will simply use L t SW , Lt , and L
S t SE , corresponding to neighbors at south west, south, and south east, respectively. We illustrate the relevant 2-time slice dynamic
Each location on the 2D grid is a random variable, representing whether that location is Empty, or occupied by a blood vessel cell. Blood vessel cells are two types: the tip of a blood vessel that has the potential to grow (henceforth called an Active Cell) or the body of the blood vessel (henceforth called the Stalk Cell). Therefore, the domain of random variable is [Active Cell, Stalk Cell, Empty], abbreviated henceforth as [AC, SC, E].
We model the conditional probability distribution, (CPD), P
| L(t+1) Lt, Lt
SW , LtS , LtSE illustrated in Figure 4.
To give a simple overview, at each step in time, an Active Cell elongates and moves into a nearby Empty location, forming the body of the blood vessel (i.e., Stalk Cell) in the process. as a tree CPD as
) t+1 • The tip of a blood vessel (AC) at time t becomes the
body (SC) at time t + 1.
P L(t+1) Lt = AC small noise parameter.
= h✏, 1 L(t+1)|Lt = AC, LtSW , LtS , LtSE = 2✏, ✏i, where ✏ is a • A
Stalk Cell at time t either continues to remain a Stalk Cell at time t + 1 or it might become Active Cell with probabilas the sprout possibility. ity to
P L(t+1) Lt = SC to sprout a new
blood vessel branch.
L(t+1)|Lt = SC, LtSW , LtS , LtSE = = h , 1 ✏, ✏i. We refer • An Empty location at time t will remain Empty at time t + 1 if none of its SW, S, or SE neighbors are Active Cell at time t; if there is an Active Cell at one or more of those neighboring locations at time t, one of them might elongate to this Empty location at time t + 1. The probability of that an Empty location being occupied by an Active Cell at time t + 1 is modeled as a Noisy-OR of its neighboring locations. iTshaatNoisisyP-ORL(to+f 1L)tS=WA, LCtS| L,LttS=E E,w,LittShWp,aLratS m,LetSteErs 0, SW , S , and SE , where 0 is leak parameter, and SW , S , and SE corresponds to the possibility that an Active Cell elongates in the NE, N, or NW direction.1 The magnitude of SW ,
S , and SE are determined by the VEGF gradient. We refer to various configurations of the parameters as the growth patterns. 4
EXPERIMENTAL SETUP,
RESULTS, AND INSIGHTS In this section, we describe the experiments we performed using various settings for the growth pattern ( ) and sprout ( ) parameters. In all the experiments to follow, we set the noise ✏ and the leak 0 parameters to 0.01. For the growth pattern, we present results for two settings: • straight-growth: h SW , S , SE i = h0.01, 0.98, 0.01i. For this pattern, the blood vessel follows a straight line, growing towards north. • uniform-growth: h SW , S , SE i = h 31 , 13 , 13 i.
For this pattern, the blood vessel has equal chance of growing towards north, north west, or north east.
For the sprout possibility, that is a Stalk Cell turning into an Active Cell, we present results for two settings: • seldom-sprout: = 0.01. For this setting, the Stalk Cell has very small chance (probability of 0.01) of becoming an Active Cell in the next time step. • always-sprout: = 0.98. For this setting, the Stalk Cell has 0.98 probability of becoming active in the next step. This is quite an unrealistic setting; we present it only for didactic purposes. We present results for four possible configurations: the cross-product of the growth patterns and sprout possibilities. We first provide detailed results on a 3 ⇥ 3 1Note that SW denotes the probability that an Active Cell at the SW of an Empty location will move to this Empty location; hence SW denotes the possibility that an Active Cell at SW moves in the NE direction to occupy an Empty location. grid over three time slices. Then, we present results on a bit larger scale, 9 ⇥ 9, over nine time slices. Finally, we present a framework where we quantify the uncertainty over the predictions on the last time slice and discuss how it is a↵ ected by the growth patterns and sprout possibilities.
For inference, in the 3⇥ 3 case, we used exact inference. For the 9 ⇥ 9 case, we used forward sampling. Note that we are able to use forward sampling in our settings because we provide the initial condition (all locations at time t = 0) as evidence and compute probabilities for the remaining time slices. 4.1
Detailed Results for 3 ⇥ 3 In this toy setting, we provide the evidence for the initial configuration of the experiment, i.e., we provide evidence for all locations for time t = 0, and compute probabilities for all locations for all future time slices. That is, we compute P (L1, L2|L0), where Lt denotes all locations at time t. For t = 0, we provide the evidence as follows: the middle of the bottom row is set as the tip of the blood vessel (i.e, L0x=1,y=0 = AC) and the rest of the locations are set as Empty. Figure 5 illustrates this setting.
E E E
E E AC
E E E The straight-growth results are presented in Figures 6 and 7, and uniform-growth results are presented in Figures 8 and 9.
The simplest setting where the blood vessel grows in a straight path and that does not sprout at all (Figure 6) is fairly straightforward to analyze. The tip of the blood vessel migrates one location towards north at each step, forming the body of the vessel along the process. This setting, therefore, serves as a sanity check. In the next setting, which is presented in Figure 7, we keep the growth pattern the same (straight-growth) but increase the sprout possibility to 0.98 (always-sprout). In this setting, the blood vessel grows towards north as expected. Unlike the seldom-sprout case, however, a Stalk Cell at time t = 1 became active at time t = 2.
Next, we present results for the uniform-growth cases. In this setting, the blood vessel has uniform AC SC AC SC .02 .01 .00 .00 .00 .01 .34 .01 .00 .00 .00 .98 .01 .00 .00 .98 .02 .01 .00 .00 .00 .01 .34 .01 .00 .00 .00 .02 .02 .01 .02 .02 .22 .02 .02 .01 .33 .02 .01 .96 .01 .96 .02 t = 2 .31 .02 .96 .01 .33 .02 t = 2 .02 .02 .01 .02 .02 .22 .02 .02 .01 .33 .02 .02 .01 .00 .00 .00 .01 .34 .01 .00 .00 .00 .98 .01 .00 .00 .98 .02 .01 .00 .00 .00 .01 .34 .01 .00 .00 .00 .02 .01 .01 .02 .02 .22 .02 .01 .01 .34 .02 .01 .01 .01 .96 .97 .02 .01 .01 .02 .02 .22 .02 .01 .01 .34 .02 probability of growing towards NW, N, and NE. In the seldom-sprout case (Figure 8), the Active Cell at t = 0 turned into a Stalk Cell at time t = 1 and remained a Stalk Cell at time t = 2. The Active Cell, unlike the straight-growth case, has equal probability of moving in all three directions. In the last time step, the middle of the top row has higher probability (.31) than the sides (.22) simply because the middle location can be reached from more locations compared to the side locations. The always-sprout case (Figure 9) is similar except a Stalk Cell at t = 1 becomes an Active Cell at t = 2.
These toy experiments provide insights into how the process typically works. Next, we present results for the 9 ⇥ 9 grid. 4.2
Summary Results for 9 ⇥ 9 Similar to the 3 ⇥ 3 grid, we provide evidence for t = 0 case and compute probabilities for the remaining eight time slices. In the initial configuration, the middle
Cell .06 .06 .06 .10 .64 .10 .06 .05 .05 .17 .17 .18 .22 .75 .22 .17 .18 .17 .06 .06 .07 .06 .01 .06 .06 .06 .06 .17 .18 .17 .18 .13 .17 .18 .18 .18 .05 .05 .05 .05 .01 .05 .05 .05 .06 .17 .16 .18 .24 .88 .24 .17 .18 .18 .05 .05 .05 .04 .01 .05 .05 .05 .05 .16 .17 .16 .16 .09 .16 .17 .17 .16 .04 .04 .04 .04 .01 .04 .04 .04 .04 .16 .14 .16 .22 .91 .22 .15 .16 .15 .04 .04 .03 .04 .01 .03 .04 .03 .04 .13 .14 .14 .13 .07 .14 .14 .14 .14 .03 .03 .03 .03 .01 .03 .03 .03 .03 .12 .11 .12 .17 .92 .16 .12 .12 .12 .02 .03 .02 .02 .01 .02 .02 .02 .02 .09 .09 .10 .10 .07 .09 .09 .10 .09 .02 .01 .01 .01 .01 .01 .01 .01 .01 .06 .06 .06 .06 .86 .05 .05 .06 .05 .22 .23 .22 .23 .23 .23 .24 .23 .23 .14 .14 .14 .15 .14 .14 .14 .14 .14 .23 .23 .23 .27 .86 .27 .23 .23 .22 .14 .14 .15 .18 .77 .18 .14 .15 .14 .23 .23 .23 .28 .87 .26 .24 .23 .22 .14 .15 .14 .14 .11 .14 .15 .14 .15 .23 .22 .23 .26 .88 .26 .23 .23 .22 .14 .13 .14 .20 .89 .20 .14 .15 .14 .21 .21 .22 .25 .89 .25 .22 .21 .21 .13 .14 .13 .13 .08 .13 .14 .14 .13 .20 .20 .20 .23 .89 .23 .20 .20 .20 .12 .11 .12 .18 .91 .18 .12 .12 .12 .18 .17 .17 .20 .90 .19 .18 .17 .18 .10 .11 .11 .10 .07 .11 .11 .10 .10 .14 .14 .14 .15 .91 .15 .14 .14 .14 .08 .08 .08 .11 .92 .10 .08 .08 .08 .10 .10 .10 .10 .93 .10 .10 .10 .10 .06 .06 .06 .06 .07 .06 .05 .06 .06 straight-growth – seldom-sprout straight-growth – always-sprout of the bottom row is set as an Active Cell and the remaining locations are set as Empty. Due to space limitations, we present results for only the last time slice, t = 8. The straight-growth case is shown in Figure 10 and the uniform-growth case is shown in Figure 11.
In the straight-growth seldom-sprout case (the left side of Figure 10), we see a straight blood vessel for the middle of the grid, where every cell of the blood vessel except the tip is a Stalk Cell and the tip is an Active Cell, as expected. In the always-sprout case (the right side of Figure 10), the Stalk Cells and Active Cells alternate, again as expected. In the uniform-growth seldom-sprout case (the left side of Figure 11), the blood vessel can be anywhere on the grid, except, as expected, the middle locations have higher probability. In the always-sprout case (the right side of Figure 11), the Stalk Cells and Active Cells alternate, as expected. Additionally, the probabilities for locations being a blood vessel (either Stalk to Active) are higher in the always-sprout case compared to the seldom-sprout case, again as expected.
The results so far have been nothing surprising, but only confirming our expectations. The value of the DBNs, however, lies at their capability to reason with spatial and temporal uncertainty as well as their potential for future directions. We discuss one of the
Cell .04 .06 .07 .08 .09 .08 .07 .06 .04 .17 .22 .24 .26 .26 .25 .24 .22 .17 .03 .04 .03 .04 .04 .04 .04 .04 .03 .16 .20 .21 .20 .20 .22 .21 .20 .16 .03 .04 .03 .04 .03 .04 .04 .03 .03 .21 .30 .34 .40 .40 .38 .34 .29 .21 .03 .04 .03 .03 .03 .03 .03 .03 .03 .15 .19 .18 .17 .17 .17 .19 .19 .15 .03 .03 .03 .03 .03 .03 .03 .03 .03 .22 .37 .51 .60 .63 .59 .50 .36 .21 .02 .03 .03 .03 .02 .03 .03 .03 .02 .15 .17 .15 .13 .12 .13 .15 .17 .15 .02 .03 .02 .02 .02 .02 .02 .03 .02 .13 .16 .52 .71 .80 .70 .53 .16 .13 .02 .02 .02 .02 .02 .02 .02 .02 .02 .10 .13 .13 .09 .09 .09 .12 .13 .10 .01 .01 .01 .01 .01 .01 .02 .02 .02 .06 .06 .06 .05 .86 .06 .06 .06 .06 .17 .19 .21 .20 .21 .21 .20 .20 .16 .13 .16 .17 .17 .16 .17 .17 .16 .13 .17 .23 .23 .26 .27 .26 .25 .23 .18 .14 .18 .21 .22 .23 .21 .21 .18 .14 .18 .22 .25 .27 .27 .27 .26 .23 .17 .13 .16 .16 .15 .17 .16 .17 .16 .13 .17 .22 .26 .28 .29 .29 .26 .22 .17 .18 .26 .34 .40 .41 .39 .33 .25 .18 .17 .22 .26 .30 .31 .30 .26 .21 .17 .12 .15 .15 .14 .13 .14 .15 .15 .12 .15 .21 .26 .32 .34 .31 .26 .20 .15 .12 .29 .47 .62 .66 .62 .46 .28 .11 .15 .16 .25 .32 .39 .33 .24 .17 .14 .11 .14 .12 .10 .09 .10 .12 .14 .11 .13 .15 .14 .40 .41 .41 .14 .14 .13 .09 .10 .10 .71 .71 .71 .10 .10 .08 .10 .10 .10 .10 .91 .10 .10 .10 .10 .06 .06 .06 .06 .07 .06 .06 .06 .05 uniform-growth – seldom-sprout uniform-growth – always-sprout future directions here supplemented with some preliminary results, and discuss more future directions in Section 6. 4.3
Quantifying Uncertainty Given an initial condition, L0, the tissue engineers are interested in the final status of the tissue, LT , where T denotes the final step of the experiment. Because real-world experiment take a long time, mostly weeks, they would like to be able to stop an experiment at time t < T and still be able to reason about time T . Therefore, they are interested in the following question: given an initial condition L0, if we stop the experiment at time t, what is the uncertainty over LT ? More practically: when is the earliest time we can stop an experiment so that the uncertainty over the last time slice is below a pre-specified target ?. It is important to note that when an experiment is stopped, the researchers dissect the tissue to analyze its properties, such as vascularization, and hence the experiment cannot continue beyond that point.
Given an uncertainty measure, this question can be formulated rather straightforwardly using DBNs. Let U N C P LT |l0, lt denote the uncertainty over the predictions over the last time slice, given the initial condition L0 = l0 and the status of the experiment at time t, Lt = lt. Then, we simply need to find argmin U N C P t<T
LT |l0, lt < Obviously, even though we know the initial conditions l0, we do not know the status of the experiment at time t > 0 unless we stop the experiment. Therefore, we need to take an expectation over all possible outcomes at time t: argmin X P t<T lt
Lt = lt|l0 U N C P
LT |l0, lt < where the subscript lt in the summation ranges over all possible configurations of Lt.
Unfortunately the number of all possible configurations for an n ⇥ n grid is 3n⇥ n, which is clearly intractable to solve. We leave a more systematic solution for future direction and present results for the case where the summation is replaced with the most probable lt|l0. For the U N C measure, there are a number of possibilities, including the entropy. We present results where we compute the conditional error of the most probable blood vessel path. That is, for the mostlikely blood vessel path, we sum 1 P (SC|lt, l0) for the body of the vessel and add 1 P (AC|lt, l0) for the tip of the blood vessel.
We experimented with the 9 ⇥ 9 grid and we set the sprout possibility to = 0.01 so that the most probably path does not have any branches. We present the uncertainty values for straight-growth and uniform-growth patterns in Figure 12. The x axis represents the time we would stop the experiment and the y axis plots the uncertainty. As expected, the uncertainty is much higher for the uniform-growth case and that uncertainty goes down for both growth patterns as we provide evidence for later time steps. We scratched only the surface of this important problem, leaving many interesting research problems for future work, some of which are discussed in Section 6. 5
RELATED
WORK
Tissue engineering experiments typically are
performed in-vivo usually on mice and in-vitro in glass
on lab. Researchers experiment with various settings
including the porosity of the sca↵ old that the tissue is
expected to hold on to, the VEGF distribution, and
initial blood vessel sprout locations [
Uncertainties for Two Growth Patterns
8.00
gineering process [
DBNs have been previously used for many
practical applications. Examples include managing
water resources [
CURRENT LIMITATIONS AND FUTURE DIRECTIONS There are two lines of work that we would like to pursue in the future. The first type is enriching the model, lifting some of the assumptions we made. The second type of work is a new line of research that we refer to as active inference, which we will describe shortly. We made a series of simplifying assumptions in our current DBN model. One such assumption is that the tissue space is 2D, whereas in reality it is obviously 3D. The 2D assumption allowed us to work with much fewer random variables. Additionally, in 2D, the number of parents for a variable is four whereas in 3D, the number of parents is ten (itself in the previous time slice and nine locations under it). It is rather straightforward to move from 2D to 3D from a representation perspective. However, scalability both in terms of memory and computational time is a challenge. Another assumption we made is that the gradient of the VEGF is fixed throughout the grid. That is, we assumed the and the values are fixed across the grid. In reality, however, the growth factor is expected to have steeper gradient when it is closer to the source of the distressed tissue cell and it is expected to be more uniform as we get further away from the distressed cell. Our simplifying assumption can be easily lifted by providing a growth factor distribution across the grid and then translating it into the necessary and parameters.
A limitation that is harder to address is scalability. In
our experiment section (Section 4), we experimented
with 3 ⇥ 3 and 9 ⇥ 9 grids. These were trivial to
experiment with. In reality, however, we need to deal with
thousands if not millions of random variables over a
much longer period of time. This will raise
scalability issues both in terms of memory and in terms of
computation time. Lifted inference [
Another line of research is to formulate and run
active inference for dynamic Bayesian networks [
CONCLUSIONS We presented a dynamic Bayesian network model for vascularization in engineered tissues. This DBN enables i) spatial and temporal reasoning for understanding of vascularization, ii) formulation and investigation of various parameter settings for vascularization, and iii) formulation of uncertainty and active information gathering to minimize uncertainty. We presented initial results that provide insights in to the vascularization process. Though the DBN model currently represents an oversimplification of the reality, it is the first and hence novel application of DBNs to vascularization. As such, it avails itself to many interesting research challenges and opportunities.
Acknowledgments This material is based upon work supported by the National Science Foundation under grant no. IIS1125412. We thank Ali Cinar, Judith Zawojewski, Eric Brey, Hamidreza Mehdizadeh, and Elif Bayrak for providing information and insights on tissue engineering.