Improved Tissue Regeneration. *Nanotoday,* Vol.4, No.1, (February 2009), pp. 66-80 **Part 4**

**Modeling and Assessment of Regeneration** 

470 Tissue Regeneration – From Basic Biology to Clinical Application

Zhang, L. & Webster, T.J. (2009). Nanotechnology and Nanomaterials: Promises for

**21** 

*Japan* 

Shogo Miyata *KEIO University* 

**Non-Invasive Evaluation Method for Cartilage** 

**Tissue Regeneration Using Quantitative-MRI** 

Articular cartilage is an avascular tissue covering articulating surfaces of bones and it functions to bear loads and reduce friction in diarthrodial joints. The cartilage can be regarded as a porous gel, mainly composed of large proteoglycan (PG) aggregates having a negative fixed-charge density (nFCD), a water-swollen network of collagen fibrils, and interstitial water, all of which play important roles in load-bearing properties (Lee et al.,

Although articular cartilage may function well over a lifetime, traumatic injury or the degenerative changes associated with osteoarthritis (OA) can significantly erode the articular layer, leading to joint pain and instability. Because of its avascular nature, articular cartilage has a very limited capacity to regenerate and repair. It is well-known that the natural response of articular cartilage to damage is variable and, at best,

Therefore, numerous studies have reported tissue-engineering approaches to restore degenerated cartilage and to repair defects; these approaches involve culturing autologous chondrocytes *in vitro* to create three-dimensional tissue that is subsequently implanted. In these tissue engineering approaches, it is important to assess the biomechanical and biochemical properties of the engineered cartilage. These material properties of the engineered constructs are detectable only via direct measurements that are invasive and require destructive treatments such as histological analysis, biochemical quantification, and mechanical testing. The application and utilization of these tissue-engineering approaches in a clinical setting requires a non-invasive method of evaluating biomechanical and biochemical properties of the actual regenerated cartilage for transplantation. Moreover, the method should be applicable to various aspects of cartilage regenerative medicine, including the characterization of the regenerated tissue during *in vitro* culture and *in vivo*

Magnetic resonance imaging (MRI) of articular cartilage is well accepted and has become common in recent years. Quantitative MRI techniques have been successfully developed to measure the macromolecular state within cartilage tissue. For example, the relationship between the water content of the degenerated cartilage and water self-diffusion has been

**1. Introduction** 

1981; Mow et al., 1980).

evaluation after transplantation.

unsatisfactory.

### **Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI**

Shogo Miyata *KEIO University Japan* 

#### **1. Introduction**

Articular cartilage is an avascular tissue covering articulating surfaces of bones and it functions to bear loads and reduce friction in diarthrodial joints. The cartilage can be regarded as a porous gel, mainly composed of large proteoglycan (PG) aggregates having a negative fixed-charge density (nFCD), a water-swollen network of collagen fibrils, and interstitial water, all of which play important roles in load-bearing properties (Lee et al., 1981; Mow et al., 1980).

Although articular cartilage may function well over a lifetime, traumatic injury or the degenerative changes associated with osteoarthritis (OA) can significantly erode the articular layer, leading to joint pain and instability. Because of its avascular nature, articular cartilage has a very limited capacity to regenerate and repair. It is well-known that the natural response of articular cartilage to damage is variable and, at best, unsatisfactory.

Therefore, numerous studies have reported tissue-engineering approaches to restore degenerated cartilage and to repair defects; these approaches involve culturing autologous chondrocytes *in vitro* to create three-dimensional tissue that is subsequently implanted. In these tissue engineering approaches, it is important to assess the biomechanical and biochemical properties of the engineered cartilage. These material properties of the engineered constructs are detectable only via direct measurements that are invasive and require destructive treatments such as histological analysis, biochemical quantification, and mechanical testing. The application and utilization of these tissue-engineering approaches in a clinical setting requires a non-invasive method of evaluating biomechanical and biochemical properties of the actual regenerated cartilage for transplantation. Moreover, the method should be applicable to various aspects of cartilage regenerative medicine, including the characterization of the regenerated tissue during *in vitro* culture and *in vivo* evaluation after transplantation.

Magnetic resonance imaging (MRI) of articular cartilage is well accepted and has become common in recent years. Quantitative MRI techniques have been successfully developed to measure the macromolecular state within cartilage tissue. For example, the relationship between the water content of the degenerated cartilage and water self-diffusion has been

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 475

assessing both biomechanical and biochemical properties during *in vitro* culture, and has been used widely in cartilage mechanobiology. Chondrocyte-seeded agarose gels were

Articular chondrocytes were obtained from the glenohumeral joints of freshly slaughtered 4 to 6-week-old calves, from a local abattoir. Articular cartilage was excised from the humeral head, diced into ~1 mm3 pieces, then shaken gently in Dulbecco's Modified Eagle's Medium/Ham's F12 (DMEM/F12) supplemented with 5% fetal bovine serum (FBS), 0.2% collagenase type II, and antibiotics-antimycotics, for 8–10 h at 37°C. Cells were then isolated from the digest by centrifugation and rinsed twice with phosphate buffered saline (PBS). Finally, after the isolated cells were resuspended with feed medium (DMEM/F12 supplemented with 20% FBS, 50 g/mL L-ascorbic acid, and antibiotics-antimycotics), and

The isolated chondrocytes in the feed medium were mixed with an equal volume of PBS containing agarose with a low melting temperature (Agarose type VII, Sigma, MO) at 37°C, to prepare 1.5 × 107 cells/mL in 2% (wt/vol) agarose gel; it was then cast in a custom-made mold to make a large gel plate. After gelling at 4°C for 25 minutes, approximately 50 disks of 8-mm diameter, 1.5-mm thickness were cored out from the large gel plate with a biopsy punch. The chondrocyte-seeded agarose disks were fed 2.5 mL feed medium/disk, every

(a) (b)

(c)

Alcian blue-stained sections of the cultured specimens are shown in Fig. 1. Over the culture period, the chondrocytes in the agarose gel appeared rounded in shape, similar to those in the "native" articular cartilage. As shown in Fig. 1, the chondrocytes synthesized a thin shell of pericellular matrix (~day 10) and expanded the volume of the cartilaginous matrix

Fig. 1. Histological appearance of tissue-engineered cartilage at day 3 (a), day 10 (b), and day 28 (c), stained with alcian blue (Miyata et al., 2010). Scale bar = 100 m.

(~day 28).

prepared as described previously (Miyata et al., 2006; Miyata et al., 2004).

the total number of cells was counted with a hemocytometer.

other day and maintained in a 5% CO2 atmosphere at 37°C.

reported (Shapiro et al., 2001), while the transverse relaxation time T2 has been related to collagen concentration (Fragonas et al., 1998) and the spatial distribution of collagen, including both fibril orientation and organization (Nieminen et al., 2001; Xia et al., 2002).

The gadolinium-diethylene triamine pentaacetic acid (Gd-DTPA2–)-enhanced T1 imaging technique has been used to predict PG content (Bashir et al., 1996) and spatial distribution (Bashir et al., 1999). Furthermore, nFCD can be estimated from consecutive T1 relaxation time measurements using Gd-DTPA2–-enhanced MRI and related to PG concentration. This MRI technique is already well-known as the "delayed Gadolinium Enhanced Magnetic Resonance Imaging of Cartilage" (dGEMRIC) technique. This technique is based on the utilization of the two-negative charge of the MRI contrast agent (i.e., Gd-DTPA2–). Sulfated glycosaminoglycans (sGAG) in the PGs are negatively charged in the cartilage, giving rise to nFCD; the electric exclusion force between this nFCD and the negatively charged contrast agent result in the inverse distribution of the contrast agent to the PG distribution in the cartilage. Consequently, relaxation time (T1) and nFCD—as determined by dGEMRIC technique —correlate with PG concentration.

Previous studies have reported that, in tissue-engineered cartilage, MR measurements of regenerated cartilage showed correlations with biochemical properties (Potter et al., 2000) and biomechanical properties (Chen et al., 2003). Additionally, the sGAG content and the compressive modulus—the latter of which was determined by unconfined compression tests—showed a trend toward correlation with the nFCD, as determined by the Gd-DTPA2–-enhanced MRI technique (Chen et al., 2003; Ramaswamy et al., 2008). In our earlier study, we reported that the nFCD of tissue-engineered cartilage determined by GD-DTPA2–-enhanced MRI has been found to correlate with sGAG content (Miyata et al., 2006).

Although the non-invasive assessment of tissue integration and the non-destructive evaluation of molecular structure of the engineered cartilage are important, we believe no previous study has fully evaluated the relationships between the biomechanical properties and MRI measurements of regenerated cartilage consisting of articular chondrocytes. Previous study has indicated that MR images of autologous chondrocyte transplants may show clinically significant variations; neither biochemical properties nor the FCD of regenerated articular cartilage has been evaluated.

In this chapter, we introduce our evaluation technique for tissue-engineered cartilage using quantitative-MRI. We tested the hypothesis that MRI measurements of tissue-engineered cartilage correlate with biomechanical and biochemical properties and that these novel approaches can be used to evaluate cartilaginous matrix material properties during tissue regeneration.

#### **2. Quantitative Magnetic Resonance Imaging (MRI) of tissue engineered cartilage**

#### **2.1 Isolation of chondrocytes and preparation of chondrocyte-seeded agarose constructs**

We used agarose gel culture for tissue-engineered cartilage model, because agarose is a biocompatible, thermosensitive hydrogel that offers superior homogeneity and stability for

reported (Shapiro et al., 2001), while the transverse relaxation time T2 has been related to collagen concentration (Fragonas et al., 1998) and the spatial distribution of collagen, including both fibril orientation and organization (Nieminen et al., 2001; Xia et al., 2002).

The gadolinium-diethylene triamine pentaacetic acid (Gd-DTPA2–)-enhanced T1 imaging technique has been used to predict PG content (Bashir et al., 1996) and spatial distribution (Bashir et al., 1999). Furthermore, nFCD can be estimated from consecutive T1 relaxation time measurements using Gd-DTPA2–-enhanced MRI and related to PG concentration. This MRI technique is already well-known as the "delayed Gadolinium Enhanced Magnetic Resonance Imaging of Cartilage" (dGEMRIC) technique. This technique is based on the utilization of the two-negative charge of the MRI contrast agent (i.e., Gd-DTPA2–). Sulfated glycosaminoglycans (sGAG) in the PGs are negatively charged in the cartilage, giving rise to nFCD; the electric exclusion force between this nFCD and the negatively charged contrast agent result in the inverse distribution of the contrast agent to the PG distribution in the cartilage. Consequently, relaxation time (T1) and nFCD—as determined by dGEMRIC

Previous studies have reported that, in tissue-engineered cartilage, MR measurements of regenerated cartilage showed correlations with biochemical properties (Potter et al., 2000) and biomechanical properties (Chen et al., 2003). Additionally, the sGAG content and the compressive modulus—the latter of which was determined by unconfined compression tests—showed a trend toward correlation with the nFCD, as determined by the Gd-DTPA2–-enhanced MRI technique (Chen et al., 2003; Ramaswamy et al., 2008). In our earlier study, we reported that the nFCD of tissue-engineered cartilage determined by GD-DTPA2–-enhanced MRI has been found to correlate with sGAG content (Miyata et al.,

Although the non-invasive assessment of tissue integration and the non-destructive evaluation of molecular structure of the engineered cartilage are important, we believe no previous study has fully evaluated the relationships between the biomechanical properties and MRI measurements of regenerated cartilage consisting of articular chondrocytes. Previous study has indicated that MR images of autologous chondrocyte transplants may show clinically significant variations; neither biochemical properties nor the FCD of

In this chapter, we introduce our evaluation technique for tissue-engineered cartilage using quantitative-MRI. We tested the hypothesis that MRI measurements of tissue-engineered cartilage correlate with biomechanical and biochemical properties and that these novel approaches can be used to evaluate cartilaginous matrix material properties during tissue

**2. Quantitative Magnetic Resonance Imaging (MRI) of tissue engineered** 

**2.1 Isolation of chondrocytes and preparation of chondrocyte-seeded agarose** 

We used agarose gel culture for tissue-engineered cartilage model, because agarose is a biocompatible, thermosensitive hydrogel that offers superior homogeneity and stability for

technique —correlate with PG concentration.

regenerated articular cartilage has been evaluated.

2006).

regeneration.

**cartilage** 

**constructs** 

assessing both biomechanical and biochemical properties during *in vitro* culture, and has been used widely in cartilage mechanobiology. Chondrocyte-seeded agarose gels were prepared as described previously (Miyata et al., 2006; Miyata et al., 2004).

Articular chondrocytes were obtained from the glenohumeral joints of freshly slaughtered 4 to 6-week-old calves, from a local abattoir. Articular cartilage was excised from the humeral head, diced into ~1 mm3 pieces, then shaken gently in Dulbecco's Modified Eagle's Medium/Ham's F12 (DMEM/F12) supplemented with 5% fetal bovine serum (FBS), 0.2% collagenase type II, and antibiotics-antimycotics, for 8–10 h at 37°C. Cells were then isolated from the digest by centrifugation and rinsed twice with phosphate buffered saline (PBS). Finally, after the isolated cells were resuspended with feed medium (DMEM/F12 supplemented with 20% FBS, 50 g/mL L-ascorbic acid, and antibiotics-antimycotics), and the total number of cells was counted with a hemocytometer.

The isolated chondrocytes in the feed medium were mixed with an equal volume of PBS containing agarose with a low melting temperature (Agarose type VII, Sigma, MO) at 37°C, to prepare 1.5 × 107 cells/mL in 2% (wt/vol) agarose gel; it was then cast in a custom-made mold to make a large gel plate. After gelling at 4°C for 25 minutes, approximately 50 disks of 8-mm diameter, 1.5-mm thickness were cored out from the large gel plate with a biopsy punch. The chondrocyte-seeded agarose disks were fed 2.5 mL feed medium/disk, every other day and maintained in a 5% CO2 atmosphere at 37°C.

Fig. 1. Histological appearance of tissue-engineered cartilage at day 3 (a), day 10 (b), and day 28 (c), stained with alcian blue (Miyata et al., 2010). Scale bar = 100 m.

Alcian blue-stained sections of the cultured specimens are shown in Fig. 1. Over the culture period, the chondrocytes in the agarose gel appeared rounded in shape, similar to those in the "native" articular cartilage. As shown in Fig. 1, the chondrocytes synthesized a thin shell of pericellular matrix (~day 10) and expanded the volume of the cartilaginous matrix (~day 28).

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 477

values of the engineered cartilage were averaged, and the results are summarized in Figure 6. T1 and Diff\* of the tissue-engineered cartilage had decreased with an increase in the culture time (Fig. 6a and 6c). On the other hand, T2 of the engineered cartilage showed considerably lower values than those of the PBS in the glass tube throughout the culture

time, and these values tended to increase slightly with the culture time (Fig. 6b).

Fig. 3. T1-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the T1 values derived from the MR images on day 1 (d), day 7 (e), and day 28

Fig. 4. T2-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the T2 values derived from the MR images on day 1 (d), day 7 (e), and day 28

(f) (Miyata et al., 2007).

(f) (Miyata et al., 2007).

#### **2.2 Magnetic Resonance Imaging (MRI) of cultured chondrocyte-seeded agarose gel**

Quantitative MRI evaluations were performed on a 2.0-T Biospec 20/30 System with a B-GA20 Gradient System (Bruker, Karlsruhe, Germany) with a maximum gradient strength of 100 mT/m. The MRI data acquisition and reconstruction were performed using the ParaVision (Bruker) software system. In all MRI experiments, three or four sheets of the disks were stacked in layers and placed into glass tubes containing phosphate buffered saline (PBS) (Fig. 2). The measured parameters included longitudinal (T1) and transverse (T2) relaxation time and water self-diffusion coefficient (Diff). A longitudinal relaxation time map (T1-map) was obtained with a short echo time (TE: 15 ms) spin-echo sequence with different repetition time values (TR: 100 ms to 15 s, 16 steps). A transverse relaxation time map (T2-map) was obtained with a long repetition time value (TR: 15 s) spin-echo sequence with different echo time values (TE: 30 ms to 450 ms, 29 steps). A diffusion coefficient map (Diff-map) was calculated from the images obtained using a conventional diffusion weighted spin-echo (SE-DWI, TR: 15 s, TE: 35 ms) sequence with different *b* values (0, 74, 275, 603, 1059 s/mm2). All sequences were performed with a field of view (FOV) of 50 × 50 mm2, matrix size 64 × 64, and slice thickness 3 mm. The values of the relaxation time (T1 and T2) and the relative diffusion coefficient (Diff\*) were calculated as the average of the specimen from the obtained T1-, T2-, and Diff-maps. The value of Diff\* (= DiffS/DiffP) was calculated by normalizing the diffusion coefficient of the sample (DiffS) by the diffusion coefficient of PBS (DiffP) around the sample. All MRI measurements were carried out with no contrast agent at room temperature (23°C).

Fig. 2. Schematic diagram of MR Imaging (Miyata et al., 2007).

Figure 3−5 shows the MRI maps of the engineered cartilage. At the first stage of the culture (day 1), T1 and Diff of the engineered cartilage showed values similar to those of the PBS around the cartilage; hence, it was difficult to distinguish the boundaries between the engineered cartilage and the bath solution (PBS) in the MRI maps (Fig. 3a and 5a). By the end of the culture (day 28), the boundaries were distinguished in both T1- and Diff-maps (Fig. 3a−3c and 5a−5c). In contrast, the boundary between the specimen and the PBS remained clear in the T2-map during the culture time (Fig. 4a−4c). The T1, T2, and Diff

**2.2 Magnetic Resonance Imaging (MRI) of cultured chondrocyte-seeded agarose gel**  Quantitative MRI evaluations were performed on a 2.0-T Biospec 20/30 System with a B-GA20 Gradient System (Bruker, Karlsruhe, Germany) with a maximum gradient strength of 100 mT/m. The MRI data acquisition and reconstruction were performed using the ParaVision (Bruker) software system. In all MRI experiments, three or four sheets of the disks were stacked in layers and placed into glass tubes containing phosphate buffered saline (PBS) (Fig. 2). The measured parameters included longitudinal (T1) and transverse (T2) relaxation time and water self-diffusion coefficient (Diff). A longitudinal relaxation time map (T1-map) was obtained with a short echo time (TE: 15 ms) spin-echo sequence with different repetition time values (TR: 100 ms to 15 s, 16 steps). A transverse relaxation time map (T2-map) was obtained with a long repetition time value (TR: 15 s) spin-echo sequence with different echo time values (TE: 30 ms to 450 ms, 29 steps). A diffusion coefficient map (Diff-map) was calculated from the images obtained using a conventional diffusion weighted spin-echo (SE-DWI, TR: 15 s, TE: 35 ms) sequence with different *b* values (0, 74, 275, 603, 1059 s/mm2). All sequences were performed with a field of view (FOV) of 50 × 50 mm2, matrix size 64 × 64, and slice thickness 3 mm. The values of the relaxation time (T1 and T2) and the relative diffusion coefficient (Diff\*) were calculated as the average of the specimen from the obtained T1-, T2-, and Diff-maps. The value of Diff\* (= DiffS/DiffP) was calculated by normalizing the diffusion coefficient of the sample (DiffS) by the diffusion coefficient of PBS (DiffP) around the sample. All MRI measurements were carried out with

no contrast agent at room temperature (23°C).

Fig. 2. Schematic diagram of MR Imaging (Miyata et al., 2007).

Figure 3−5 shows the MRI maps of the engineered cartilage. At the first stage of the culture (day 1), T1 and Diff of the engineered cartilage showed values similar to those of the PBS around the cartilage; hence, it was difficult to distinguish the boundaries between the engineered cartilage and the bath solution (PBS) in the MRI maps (Fig. 3a and 5a). By the end of the culture (day 28), the boundaries were distinguished in both T1- and Diff-maps (Fig. 3a−3c and 5a−5c). In contrast, the boundary between the specimen and the PBS remained clear in the T2-map during the culture time (Fig. 4a−4c). The T1, T2, and Diff values of the engineered cartilage were averaged, and the results are summarized in Figure 6. T1 and Diff\* of the tissue-engineered cartilage had decreased with an increase in the culture time (Fig. 6a and 6c). On the other hand, T2 of the engineered cartilage showed considerably lower values than those of the PBS in the glass tube throughout the culture time, and these values tended to increase slightly with the culture time (Fig. 6b).

Fig. 3. T1-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the T1 values derived from the MR images on day 1 (d), day 7 (e), and day 28 (f) (Miyata et al., 2007).

Fig. 4. T2-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the T2 values derived from the MR images on day 1 (d), day 7 (e), and day 28 (f) (Miyata et al., 2007).

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 479

(a)

(b)

(c)

Fig. 6. Longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion coefficient (c) of the tissue-engineered cartilage during the culture time (Miyata et al., 2007).

For 'native' articular cartilage, the gadolinium-diethylene triamine pentaacetic acid (Gd-DTPA2-) -enhanced T1 imaging technique has been used to predict the PG content (Bashir et

**2.3 Evaluation of fixed charge density of tissue-engineered cartilage** 

The values represent mean +/– S.D. (n = 3).

Consistent with the results of previous studies (Chen et al., 2003; Potter et al., 1998), our results showed that both T1 and Diff\* decreased with an increase in the culture time. On the other hand, T2 tended to increase slightly by the end of the culture time. To understand the MRI properties of the tissue water protons, we have to understand the behavior of water molecules in the tissue at different stages of tissue maturity. With tissue growth and development, proteoglycan and collagen molecules accumulate in the agarose gel, resulting in a large fraction of macromolecule-associated water, which is known as "bound" water. Generally, the water molecules in the "bound" condition show short T1 and T2 relaxation times due to a reduced mobility as compared to "free" water. Thus, water proton relaxation curves, which were described by a single exponential, are derived from the weighted sum of the relaxation behavior of the "free" and "bound" water molecules in the engineered cartilage. This is consistent with our results that the T1 relaxation time and Diff\* decreased with an increase in the content of cartilaginous matrix in the agarose gel. In the case of transverse relaxation, the T2 relaxation time of the engineered cartilage showed a value similar to that of the "native" articular cartilage (75–90 ms measured by our MRI system) from the early phase of the culture; further, T2 tended to increase slightly with tissue maturation. Based on this result, we speculate that the transverse relaxation of the water molecules in the engineered construct might be mainly affected by its association with the agarose molecules.

Fig. 5. Diff-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the Diff values derived from the MR images on day 1 (d), day 7 (e), and day 28 (f) (Miyata et al., 2007).

Consistent with the results of previous studies (Chen et al., 2003; Potter et al., 1998), our results showed that both T1 and Diff\* decreased with an increase in the culture time. On the other hand, T2 tended to increase slightly by the end of the culture time. To understand the MRI properties of the tissue water protons, we have to understand the behavior of water molecules in the tissue at different stages of tissue maturity. With tissue growth and development, proteoglycan and collagen molecules accumulate in the agarose gel, resulting in a large fraction of macromolecule-associated water, which is known as "bound" water. Generally, the water molecules in the "bound" condition show short T1 and T2 relaxation times due to a reduced mobility as compared to "free" water. Thus, water proton relaxation curves, which were described by a single exponential, are derived from the weighted sum of the relaxation behavior of the "free" and "bound" water molecules in the engineered cartilage. This is consistent with our results that the T1 relaxation time and Diff\* decreased with an increase in the content of cartilaginous matrix in the agarose gel. In the case of transverse relaxation, the T2 relaxation time of the engineered cartilage showed a value similar to that of the "native" articular cartilage (75–90 ms measured by our MRI system) from the early phase of the culture; further, T2 tended to increase slightly with tissue maturation. Based on this result, we speculate that the transverse relaxation of the water molecules in the engineered construct might be mainly affected by its association with the

Fig. 5. Diff-maps of day 1 (a), day 7 (b), and day 28 (c) post-inoculation specimens, and histograms of the Diff values derived from the MR images on day 1 (d), day 7 (e), and day

agarose molecules.

28 (f) (Miyata et al., 2007).

Fig. 6. Longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion coefficient (c) of the tissue-engineered cartilage during the culture time (Miyata et al., 2007). The values represent mean +/– S.D. (n = 3).

#### **2.3 Evaluation of fixed charge density of tissue-engineered cartilage**

For 'native' articular cartilage, the gadolinium-diethylene triamine pentaacetic acid (Gd-DTPA2-) -enhanced T1 imaging technique has been used to predict the PG content (Bashir et

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 481

Fig. 7. Schematic diagram of gadolinium-enhanced MRI. In all MRI measurements, the cultured specimens were put into glass tubes filled with phosphate buffered saline (PBS) or

(a) (b) (c) Fig. 8. Quantitative water proton T1 maps in the presence of Gd-DTPA2– at day 3 (a), day 7

Fig. 9. Tissue fixed-charge density, with time in culture, for tissue-engineered cartilage

(Miyata et al., 2010). \* indicates significant difference from day 0 (P < 0.05).

1 mM Gd DTPA2– (Miyata et al., 2010).

(b), and day 28 (c) (Miyata et al., 2010).

al., 1996) and spatial distribution (Bashir et al., 1999). Furthermore, the negative fixed charge density (nFCD) can be estimated from consecutive T1 relaxation time measurement using Gd-DTPA2--enhanced MRI and be related to the PG concentration. In this study, we used this dGEMRIC technique to monitor and evaluate tissue integration of the engineered cartilage.

The MRI measurements were performed with a 2.0-Tesla Bruker Biospec 20/30 system using Gd-DTPA2- contrast agent. In all MRI measurements, the specimens were put into glass tubes filled with PBS (Fig. 7). The longitudinal relaxation time map, T1-map, was obtained with a short-echo time (TE: 15 ms), spin-echo sequence with different repetition time values (TR: 100 ms to 15 s, 16 steps). Subsequently, the specimens were balanced in PBS containing 1 mM Gd-DTPA2– (Magnevist®, Nihon Schering, Osaka, Japan) for 10−12 hours; the longitudinal relaxation time map in the contrast agent, T1Gd-map, was obtained again with a short-echo time (TE: 15 ms), spin-echo sequence with different repetition time values (TR: 30 ms to 5 s, 13 steps). Finally, using the relaxivity (R) value of Gd-DTPA2– in saline (5.24 in our MRI system), the concentration of the contrast agent was estimated using the formula [Gd-DTPA2–] = 1/R(1/T1Gd – 1/T1). The negative fixed charge density (FCD) was calculated as follows

$$\text{nFCD} = \frac{[\text{Na}^+]\_b \sqrt{[\text{Gd} - \text{DTPA}^2 \text{ }^-]\_t}}{\sqrt{[\text{Gd} - \text{DTPA}^2 \text{ }^-]\_b}} - \frac{[\text{Na}^+]\_b \sqrt{[\text{Gd} - \text{DTPA}^2 \text{ }^-]\_b}}{\sqrt{[\text{Gd} - \text{DTPA}^2 \text{ }^-]\_t}} \tag{1}$$

where subscript *b* stands for bath solution and subscript *t* stands for cartilaginous tissue(Bashir et al., 1996). All MRI measurements were performed at room temperature 23°C.

In the gadolinium-enhanced MR imaging measurements, longitudinal relaxation time of the bulk PBS containing Gd-DTPA reagent showed 0.179 ± 0.06 seconds in our MRI system. The T1Gd of the cultured specimen increased as a function of tissue maturation (0.197 ± 0.001 to 0.222 ± 0.003 seconds). At the first stage of the culture (day 3), T1Gd of the tissue-engineered cartilage showed values proximate to those of the PBS containing the Gd-DTPA2– agent around the engineered cartilage; hence, it was difficult to distinguish the boundaries between the engineered cartilage and the bath solution in the T1Gd-maps (Fig. 8a). By the end of the culture (day 28), the boundaries had become distinct in the T1Gd-maps (Fig. 8). The [Gd-DTPA2–] in the engineered cartilage decreased with increases in culture time. The nFCD, as determined from the [Gd-DTPA2–] in the specimen and bath solution, increased with culture time (Fig. 9).

As time in culture lengthened, the gross appearance of the cultured disk became increasingly opaque. The DMMB assay (Farndale et al., 1986) revealed that the sGAG content of the chondrocyte/agarose disks increased as a function of tissue maturation (0.19 ± 0.27 to 13.2 ± 1.9 mg/mL-disk-vol). Finally, the sGAG content of the reconstructed cartilaginous disk reached approximately 20% of the "native" articular cartilage (data not shown).

To correlate gadolinium-enhanced MRI and biochemical properties, the sGAG content of the tissue was plotted as a function of the FCD. From the linear regression analysis, the FCD correlated significantly with the sGAG content (r = 0.95, n = 30, P < 0.001) (Fig. 10), and the tissue [Gd-DTPA2–] correlated with the sGAG content by r = 0.83, n = 30, P < 0.001.

al., 1996) and spatial distribution (Bashir et al., 1999). Furthermore, the negative fixed charge density (nFCD) can be estimated from consecutive T1 relaxation time measurement using Gd-DTPA2--enhanced MRI and be related to the PG concentration. In this study, we used this dGEMRIC technique to monitor and evaluate tissue integration of the engineered cartilage.

The MRI measurements were performed with a 2.0-Tesla Bruker Biospec 20/30 system using Gd-DTPA2- contrast agent. In all MRI measurements, the specimens were put into glass tubes filled with PBS (Fig. 7). The longitudinal relaxation time map, T1-map, was obtained with a short-echo time (TE: 15 ms), spin-echo sequence with different repetition time values (TR: 100 ms to 15 s, 16 steps). Subsequently, the specimens were balanced in PBS containing 1 mM Gd-DTPA2– (Magnevist®, Nihon Schering, Osaka, Japan) for 10−12 hours; the longitudinal relaxation time map in the contrast agent, T1Gd-map, was obtained again with a short-echo time (TE: 15 ms), spin-echo sequence with different repetition time values (TR: 30 ms to 5 s, 13 steps). Finally, using the relaxivity (R) value of Gd-DTPA2– in saline (5.24 in our MRI system), the concentration of the contrast agent was estimated using the formula [Gd-DTPA2–] = 1/R(1/T1Gd – 1/T1). The negative fixed charge density (FCD) was

2 2 [Na ] [Gd DTPA ] [Na ] [Gd DTPA ]

where subscript *b* stands for bath solution and subscript *t* stands for cartilaginous tissue(Bashir et al., 1996). All MRI measurements were performed at room temperature

In the gadolinium-enhanced MR imaging measurements, longitudinal relaxation time of the bulk PBS containing Gd-DTPA reagent showed 0.179 ± 0.06 seconds in our MRI system. The T1Gd of the cultured specimen increased as a function of tissue maturation (0.197 ± 0.001 to 0.222 ± 0.003 seconds). At the first stage of the culture (day 3), T1Gd of the tissue-engineered cartilage showed values proximate to those of the PBS containing the Gd-DTPA2– agent around the engineered cartilage; hence, it was difficult to distinguish the boundaries between the engineered cartilage and the bath solution in the T1Gd-maps (Fig. 8a). By the end of the culture (day 28), the boundaries had become distinct in the T1Gd-maps (Fig. 8). The [Gd-DTPA2–] in the engineered cartilage decreased with increases in culture time. The nFCD, as determined from the [Gd-DTPA2–] in the specimen and bath solution, increased

As time in culture lengthened, the gross appearance of the cultured disk became increasingly opaque. The DMMB assay (Farndale et al., 1986) revealed that the sGAG content of the chondrocyte/agarose disks increased as a function of tissue maturation (0.19 ± 0.27 to 13.2 ± 1.9 mg/mL-disk-vol). Finally, the sGAG content of the reconstructed cartilaginous disk reached approximately 20% of the "native" articular cartilage (data not

To correlate gadolinium-enhanced MRI and biochemical properties, the sGAG content of the tissue was plotted as a function of the FCD. From the linear regression analysis, the FCD correlated significantly with the sGAG content (r = 0.95, n = 30, P < 0.001) (Fig. 10), and the

tissue [Gd-DTPA2–] correlated with the sGAG content by r = 0.83, n = 30, P < 0.001.

2 2 [Gd DTPA ] [Gd DTPA ] *b tb b*

*b t*

(1)

calculated as follows

with culture time (Fig. 9).

23°C.

shown).

nFCD

Fig. 7. Schematic diagram of gadolinium-enhanced MRI. In all MRI measurements, the cultured specimens were put into glass tubes filled with phosphate buffered saline (PBS) or 1 mM Gd DTPA2– (Miyata et al., 2010).

Fig. 8. Quantitative water proton T1 maps in the presence of Gd-DTPA2– at day 3 (a), day 7 (b), and day 28 (c) (Miyata et al., 2010).

Fig. 9. Tissue fixed-charge density, with time in culture, for tissue-engineered cartilage (Miyata et al., 2010). \* indicates significant difference from day 0 (P < 0.05).

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 483

(a)

(b) Fig. 11. Equilibrium compressive modulus *E*eq (a) and dynamic compressive modulus *E*dyn (b), with time in culture, for the cultured chondrocyte/agarose disks (Miyata et al., 2010).

To determine the correlations between the quantitative MRI measurements and the biomechanical and biochemical properties of the tissue-engineered cartilage, we performed linear regression analyses among the MRI-derived parameters (T1, T2, Diff, and FCD), the biochemical composition (sGAG content), and the biomechanical properties (*E*eq, *E*dyn) of the

To confirm the correlation, the *E*eq of the engineered cartilage were plotted as functions of the T1, T2, and Diff, respectively. The *E*eq of the engineered cartilage (Fig. 12a) showed a strong correlation with T1 and Diff but a weak correlation with T2 (Fig. 12b and 12c). Similarly, the tissue sGAG concentration (Fig. 13a and 13c) and were found to be strongly correlated with T1 and Diff. Consistent with the results of the previous investigation (Potter et al., 2000), our results showed that T1 relaxation time and Diff showed a significant correlation with the biomechanical properties and the sGAG content of the tissueengineered cartilage. The results of recent studies have shown that the articular cartilage

**2.5 Relationships between MRI measurements and biomechanical properties of** 

\* indicates significant difference from day 0 (P < 0.05).

**cultured chondrocyte-seeded constructs** 

engineered cartilage.

Fig. 10. Scatter plots relating the tissue fixed charge density (FCD) to the sulfated glycosaminoglycan (sGAG) content (Miyata et al., 2006).

#### **2.4 Static and dynamic biomechanical testing of cultured agarose/chondrocyte constructs**

Mechanical testing of the disk-shaped specimens was performed with unconfined compression, using impermeable stainless platens in PBS at room temperature. Static compressive properties were measured in a custom-made chamber attached to a material testing device (Autograph 5kNG, Shimadzu, Kyoto, Japan). Stress relaxation tests were performed by applying a ramp displacement at 0.05 mm/min to a 20% static compressive strain, followed by relaxation to equilibrium (2,400 s). The equilibrium compressive modulus (*E*eq) was calculated from the imposed compressive strain and the equilibrium load, divided by the cross-sectional area of the specimen.

Dynamic compression tests were carried out using a viscoelastic spectrometer (DDV-MF, A&D, Tokyo, Japan) (Miyata et al., 2005). For preconditioning, a 20% static compressive strain was loaded and a sinusoidal displacement of 0.5% compressive strain was then superimposed at a frequency of 1 Hz. After equilibrium had been reached (approximately 20 min), a sinusoidal displacement of 0.5% compressive strain was applied at frequencies ranging from 0.01 to 5.0 Hz. The dynamic compressive modulus (*E*dyn) was calculated from the ratio of the measured stress amplitude and the applied strain amplitude.

Figure 11 shows the means and standard deviations of the equilibrium compressive modulus *E*eq and dynamic compressive modulus *E*dyn versus time in culture, for the tissueengineered cartilage. With respect to the static compressive property, significant differences were observed in the equilibrium compressive modulus. With increases in culture time, the *E*eq of the specimens increased and reached approximately 10% of that of "native" cartilage from which the chondrocytes were harvested (0.45 ± 0.12 MPa, n =3). With respect to the dynamic compressive property, significant differences were also observed for testing conditions. The dynamic compressive modulus *E*dyn of the engineered cartilage depended on both testing frequency and culture time. For each time point, *E*dyn increased nonlinearly with increases in frequency. The engineered cartilage also exhibited marked stiffening with time in culture. The value of *E*dyn increased with culture time at each testing frequency (0.01– 2.0 Hz).

Fig. 10. Scatter plots relating the tissue fixed charge density (FCD) to the sulfated

**2.4 Static and dynamic biomechanical testing of cultured agarose/chondrocyte** 

Mechanical testing of the disk-shaped specimens was performed with unconfined compression, using impermeable stainless platens in PBS at room temperature. Static compressive properties were measured in a custom-made chamber attached to a material testing device (Autograph 5kNG, Shimadzu, Kyoto, Japan). Stress relaxation tests were performed by applying a ramp displacement at 0.05 mm/min to a 20% static compressive strain, followed by relaxation to equilibrium (2,400 s). The equilibrium compressive modulus (*E*eq) was calculated from the imposed compressive strain and the equilibrium

Dynamic compression tests were carried out using a viscoelastic spectrometer (DDV-MF, A&D, Tokyo, Japan) (Miyata et al., 2005). For preconditioning, a 20% static compressive strain was loaded and a sinusoidal displacement of 0.5% compressive strain was then superimposed at a frequency of 1 Hz. After equilibrium had been reached (approximately 20 min), a sinusoidal displacement of 0.5% compressive strain was applied at frequencies ranging from 0.01 to 5.0 Hz. The dynamic compressive modulus (*E*dyn) was calculated from

Figure 11 shows the means and standard deviations of the equilibrium compressive modulus *E*eq and dynamic compressive modulus *E*dyn versus time in culture, for the tissueengineered cartilage. With respect to the static compressive property, significant differences were observed in the equilibrium compressive modulus. With increases in culture time, the *E*eq of the specimens increased and reached approximately 10% of that of "native" cartilage from which the chondrocytes were harvested (0.45 ± 0.12 MPa, n =3). With respect to the dynamic compressive property, significant differences were also observed for testing conditions. The dynamic compressive modulus *E*dyn of the engineered cartilage depended on both testing frequency and culture time. For each time point, *E*dyn increased nonlinearly with increases in frequency. The engineered cartilage also exhibited marked stiffening with time in culture. The value of *E*dyn increased with culture time at each testing frequency (0.01–

the ratio of the measured stress amplitude and the applied strain amplitude.

glycosaminoglycan (sGAG) content (Miyata et al., 2006).

load, divided by the cross-sectional area of the specimen.

**constructs** 

2.0 Hz).

Fig. 11. Equilibrium compressive modulus *E*eq (a) and dynamic compressive modulus *E*dyn (b), with time in culture, for the cultured chondrocyte/agarose disks (Miyata et al., 2010). \* indicates significant difference from day 0 (P < 0.05).

#### **2.5 Relationships between MRI measurements and biomechanical properties of cultured chondrocyte-seeded constructs**

To determine the correlations between the quantitative MRI measurements and the biomechanical and biochemical properties of the tissue-engineered cartilage, we performed linear regression analyses among the MRI-derived parameters (T1, T2, Diff, and FCD), the biochemical composition (sGAG content), and the biomechanical properties (*E*eq, *E*dyn) of the engineered cartilage.

To confirm the correlation, the *E*eq of the engineered cartilage were plotted as functions of the T1, T2, and Diff, respectively. The *E*eq of the engineered cartilage (Fig. 12a) showed a strong correlation with T1 and Diff but a weak correlation with T2 (Fig. 12b and 12c). Similarly, the tissue sGAG concentration (Fig. 13a and 13c) and were found to be strongly correlated with T1 and Diff. Consistent with the results of the previous investigation (Potter et al., 2000), our results showed that T1 relaxation time and Diff showed a significant correlation with the biomechanical properties and the sGAG content of the tissueengineered cartilage. The results of recent studies have shown that the articular cartilage

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 485

(a)

(b)

(c) Fig. 13. Scatter plots for the relationship between the equilibrium compressive modulus *E*eq and longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion

One possible explanation is that the changes in the biophysical properties might be mainly due to the altered sGAG content, and the synthesis of collagen and the reorganization of

coefficient (c) (Miyata et al., 2007). Solid line represents the linear regression line.

collagen network might be insufficient in the agarose gel culture.

degeneration induced by collagenase treatment resulted in changes in T2 relaxation time and the equilibrium modulus (Nieminen et al., 2000). In the present study, slight increase in T2 values was observed during the culture.

Fig. 12. Scatter plots for the relationship between the equilibrium compressive modulus *E*eq and longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion coefficient (c) (Miyata et al., 2007). Solid line represents the linear regression line.

degeneration induced by collagenase treatment resulted in changes in T2 relaxation time and the equilibrium modulus (Nieminen et al., 2000). In the present study, slight increase in

(a)

(b)

(c) Fig. 12. Scatter plots for the relationship between the equilibrium compressive modulus *E*eq and longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion

coefficient (c) (Miyata et al., 2007). Solid line represents the linear regression line.

T2 values was observed during the culture.

Fig. 13. Scatter plots for the relationship between the equilibrium compressive modulus *E*eq and longitudinal relaxation time (a), transverse relaxation time (b), and relative diffusion coefficient (c) (Miyata et al., 2007). Solid line represents the linear regression line.

One possible explanation is that the changes in the biophysical properties might be mainly due to the altered sGAG content, and the synthesis of collagen and the reorganization of collagen network might be insufficient in the agarose gel culture.

Non-Invasive Evaluation Method for Cartilage Tissue Regeneration Using Quantitative-MRI 487

crucial role in static compressive behavior, while collagen bears a dynamic compressive load (Korhonen et al., 2003). Therefore, the nFCD—which is to say, the sGAG content—might show a higher correlation with the equilibrium modulus than with the dynamic modulus. Moreover, the dynamic modulus showed a trend toward correlation with the nFCD at lower frequencies than that of higher frequencies. That might reflect the collagen network levels regenerated in the agarose gel. Nonetheless, the results of recent studies have shown that variations in collagen architecture among varieties of articular cartilage decreased the significance of correlations between Gd-DTPA2–-enhanced MRI and mechanical properties, because the architecture of the collagen network, as well as PGs, plays an important role in the mechanicalproperties of articular cartilage (Nissi et al., 2007). In the present study, the chondrocytes in agarose gel reconstructed the immature collagen network, prompting a low-level effect on the compressive property compared to "native" articular cartilage; therefore, significant correlations might be found between Gd-DTPA2–-enhanced MRI and compressive properties. From these facts, our evaluation methods using Gd-DTPA2–-

In conclusion, we evaluated the changes in the quantitative MRI parameters and matrix FCD of tissue-engineered cartilage that consisted of articular chondrocytes and hydrogels. We found significant linear correlations between the quantitative MRI measurements and the biomechanical and biochemical properties of the engineered cartilage. Finally, we suggest that the quantitative MRI technique can be a useful, non-invasive approach to evaluate the

This research was supported in part by a Grant-in-Aid for Young Scientists (B) (No.

Bashir, A.; Gray, M. L. & Burstein, D. (1996). Gd-DTPA2- as a Measure of Cartilage

Bashir, A.; Gray, M.L.; Hartke, J. & Burstein, D. (1999). Nondestructive Imaging of Human

Chen, C. T.; Fishbein, K.W.; Torzilli, P. A.; Hilger, A.; Spencer, R. G. & Horton, W. E., Jr.

Fragonas, E.; Mlynárik, V.; Jellús, V.; Micali, F.; Piras, A.; Toffanin, R.; Rizzo, R. & Vittur, F.

Cartilage Glycosaminoglycan Concentration by MRI. *Magnetic Resonance in* 

(2003). Matrix Fixed-Charge Density as Determined by Magnetic Resonance Microscopy of Bioreactor-Derived Hyaline Cartilage Correlates with Biochemical and Biomechanical Properties. *Arthritis and Rheumatism*, Vol.48, pp.1047-1056 Farndale, R. W.; Buttle, D. J. & Barrett, A. J. (1986). Improved Quantitation and

Discrimination of Sulphated Glycosaminoglycans by Use of Dimethylmethylene

(1998). Correlation between Biochemical Composition and Magnetic Resonance Appearance of Articular Cartilage. *Osteoarthritis and Cartilage*, Vol.6, pp.24–32 Korhonen, R.K.; Laasanen, M.S.; Töyräs, J.; Lappalainen, R.; Helminen, H.J. & Jurvelin, J.S.

(2003). Fibril Reinforced Poroelastic Model Predicts Specifically Mechanical

enhanced MRI could be applicable at the earlier stage of tissue regeneration.

biomechanical properties of regenerated cartilage during *in vitro* culturing process.

18700414) from the Ministry of Education, Science, Sports and Culture of Japan.

Degradation. *Magnetic Resonance Medicine*, Vol.36, pp.665-673

Blue. *Biochimica et Biophysica Acta*, Vol.883, pp.173-177

**3. Conclusion** 

**4. Acknowledgment** 

*Medicine*, Vol.41, pp.857–865

**5. References** 

Fig. 14. Typical scatter plots relating the tissue fixed-charge density to equilibrium compressive modulus *E*eq (a) and dynamic compressive modulus *E*dyn at 0.5 Hz (b) (Miyata et al., 2010).


Table 1. Linear Pearson correlations between biomechanical and Gd-DTPA2–-enhanced MRI parameters in tissue-engineered cartilage (Miyata et al., 2010).

To evaluate the relationship between Gd-DTPA2–-enhanced MRI parameters and biomechanical properties, the *E*eq and *E*dyn of the engineered cartilage were plotted as functions of the nFCD, respectively. From the linear Pearson correlation analysis, it was found that nFCD correlated significantly with *E*eq and *E*dyn (Table 1, Fig. 14). The equilibrium compressive modulus showed a higher correlation than the dynamic compressive modulus of all testing frequencies, and the dynamic compressive modulus tended to show a slightly higher correlation at low frequencies (0.01–0.05 Hz). The sGAG of articular cartilage plays a crucial role in static compressive behavior, while collagen bears a dynamic compressive load (Korhonen et al., 2003). Therefore, the nFCD—which is to say, the sGAG content—might show a higher correlation with the equilibrium modulus than with the dynamic modulus. Moreover, the dynamic modulus showed a trend toward correlation with the nFCD at lower frequencies than that of higher frequencies. That might reflect the collagen network levels regenerated in the agarose gel. Nonetheless, the results of recent studies have shown that variations in collagen architecture among varieties of articular cartilage decreased the significance of correlations between Gd-DTPA2–-enhanced MRI and mechanical properties, because the architecture of the collagen network, as well as PGs, plays an important role in the mechanicalproperties of articular cartilage (Nissi et al., 2007). In the present study, the chondrocytes in agarose gel reconstructed the immature collagen network, prompting a low-level effect on the compressive property compared to "native" articular cartilage; therefore, significant correlations might be found between Gd-DTPA2–-enhanced MRI and compressive properties. From these facts, our evaluation methods using Gd-DTPA2– enhanced MRI could be applicable at the earlier stage of tissue regeneration.

#### **3. Conclusion**

486 Tissue Regeneration – From Basic Biology to Clinical Application

(a)

(b)

compressive modulus *E*eq (a) and dynamic compressive modulus *E*dyn at 0.5 Hz (b) (Miyata

Table 1. Linear Pearson correlations between biomechanical and Gd-DTPA2–-enhanced MRI

To evaluate the relationship between Gd-DTPA2–-enhanced MRI parameters and biomechanical properties, the *E*eq and *E*dyn of the engineered cartilage were plotted as functions of the nFCD, respectively. From the linear Pearson correlation analysis, it was found that nFCD correlated significantly with *E*eq and *E*dyn (Table 1, Fig. 14). The equilibrium compressive modulus showed a higher correlation than the dynamic compressive modulus of all testing frequencies, and the dynamic compressive modulus tended to show a slightly higher correlation at low frequencies (0.01–0.05 Hz). The sGAG of articular cartilage plays a

Fig. 14. Typical scatter plots relating the tissue fixed-charge density to equilibrium

 R2 *P*  FCD vs. *E*eq 0.81 < 0.001 FCD vs. *E*dyn, 0.01 Hz 0.79 < 0.001 FCD vs. *E*dyn, 0.02 Hz 0.73 < 0.001 FCD vs. *E*dyn, 0.05 Hz 0.73 < 0.001 FCD vs. *E*dyn, 0.5 Hz 0.70 < 0.001 FCD vs. *E*dyn, 2.0 Hz 0.71 < 0.001

parameters in tissue-engineered cartilage (Miyata et al., 2010).

et al., 2010).

In conclusion, we evaluated the changes in the quantitative MRI parameters and matrix FCD of tissue-engineered cartilage that consisted of articular chondrocytes and hydrogels. We found significant linear correlations between the quantitative MRI measurements and the biomechanical and biochemical properties of the engineered cartilage. Finally, we suggest that the quantitative MRI technique can be a useful, non-invasive approach to evaluate the biomechanical properties of regenerated cartilage during *in vitro* culturing process.

#### **4. Acknowledgment**

This research was supported in part by a Grant-in-Aid for Young Scientists (B) (No. 18700414) from the Ministry of Education, Science, Sports and Culture of Japan.

#### **5. References**


**0**

**22**

*The Netherlands*

**A Mathematical Model for Wound Contraction and Angiogenesis**

*Delft Institute of Applied Mathematics, Delft University of Technology*

Cutaneous wounds, ulcers or burns, as a result of external damage, such as intensive solar exposure or accidents, occur in high numbers and may cause ever lasting traumas. In some cases, the wounds are very painful and impair an individual's daily life in terms of health, and social interaction with his or her surroundings, but also in terms of mobility and ability to work or to perform leisure activities. Furthermore, the aesthetic feeling of patients is highly improved by the treatment of burns and other wounds since otherwise the wounds can make a detrimental impact on the appearance of the patient's skin. Besides aesthetic appearance, also the functioning of the damaged skin is different since infections are possibly more likely to occur near the (scar tissue of the) burn. Furthermore, the conception of temperature and pain by a patient is also altered by the burn, ulcer or cutaneous wound. Therefore, the treatment of these defects, which occurs in very high numbers, is important and hence it is crucial to

Nowadays many different treatments, such as, in more serious cases, skin transplantation, injection of cell cultures and the implant of artificial skin, are commonly carried out by plastical surgeons in order to relieve the patient and to minimize the aesthetic effects of scar tissue as much as possible. In order to be able to improve medical treatment of scars, it is important to understand the biological mechanism behind scar formation in terms of plastic deformation and a regenerated excess of collageneous matrix. More information about the treatment of scars and the biology behind the scar formation can, among others, be found in Bloemen (2011); van der Veer (2009). The deformations are caused by the contractile mechanism due to the pulling behavior of the (myo-)fibroblasts. Dallon (2008) experimentally found that the nature of contraction is predominantly determined by the following two



These mechanisms were simulated by Dallon (2010) by the use of a spring model. Since Dallon's two-dimensional model tracks the forces exerted by each individual philapod or

improve treatments and to make treatments more efficient.

**1. Introduction**

mechanisms:

matrix.

onto the collageneous fibres;

Fred Vermolen and Olmer van Rijn

Behavior of Normal, Proteoglycan Depleted and Collagen Degraded Articular Cartilage. *Journal of Biomechanics*, Vol.36, pp.1373–1379


### **A Mathematical Model for Wound Contraction and Angiogenesis**

Fred Vermolen and Olmer van Rijn *Delft Institute of Applied Mathematics, Delft University of Technology The Netherlands*

#### **1. Introduction**

488 Tissue Regeneration – From Basic Biology to Clinical Application

Lee, R. C.; Frank, E. H.; Grodzinsky, A. J. & Roylance, D. K. (1981). Oscillatory Compressional

Miyata, S.; Tateishi, T.; Furukawa, K. & Ushida, T. (2005). Influence of Structure and

Miyata, S.; Numano, T.; Homma, H.; Tateishi, T. & Ushida, T. (2007). Feasibility of

Mow, V.C.; Kuei, S. C.; Lai, W.M. & Armstrong, C.G. (1980). Biphasic Creep and Stress

Nieminen, M.T.; Rieppo, J.; Toyras, J.; Hakumaki, J.M.; Silvennoinen, J.; Hyttinen, M.M.;

Light Microscopic Study. *Magnetic Resonance in Medicine*, Vol.46, pp.487-493 Nieminen, M.T.; Toyras, J.; Rieppo, J.; Hakumaki, J.M.; Silvennoinen, J.; Helminen, H.J. &

Potter, K.; Butler, J.J.; Adams, C.; Fishbein, K.W.; McFarland, E.W.; Horton, W.E. & Spencer,

Ramaswamy, S.; Uluer, M.C.; Leen, S.; Bajaj, P.; Fishbein, K.W. & Spencer, R.G. (2008).

Shapiro, E.M.; Borthakur, A.; Kaufman, J.H.; Leigh, J.S. & Reddy, R. (2001). Water

Magnetic Resonance Imaging. *Osteoarthritis Cartilage*, Vol.9, pp.533-538 Xia, Y.; Moody, J.B. & Alhadlaq, H. (2002). Orientational Dependence of T2 Relaxation in

Magnetic Resonance Microscopy. *Matrix Biology*, Vol.17, pp.513-523 Potter, K.; Butler, J.J.; Horton, W.E. & Spencer, R.G. (2000). Response of Engineered Cartilage

Articular Cartilage. *Magnetic Resonance in Medicine*, Vol.43, pp.676-681. Nissi, M.J.; Rieppo, J.; Töyräs, J.; Laasanen, M.S.; Kiviranta, I.; Nieminen, M.T. & Jurvelin,

Maturation. *Osteoarthritis and Cartilage*, Vol.15, pp.24–32

Microscopy. *Arthritis and Rheumatism*, Vol.43, pp.1580-1590

*Engineering Part C: Methods*, Vol.14, pp.243–249

*Medicine*, Vol.48, pp.460–469

by Using Quantitative MRI. *Journal of Biomechnaics*, Vol.40, pp. 2990-2998 Miyata, S.; Homma, H.; Numano, T.; Tateishi, T. & Ushida, T. (2010). Evaluation of Negative

Mechanical Function. *JSME International Journal: C*, Vol.48, pp.547–554 Miyata, S.; Homma, H.; Numano, T.; Furukawa, K; Tateishi, T. & Ushida, T. (2006).

Cartilage. *Journal of Biomechanics*, Vol.36, pp.1373–1379

*Journal of Biomechanical Engineering*, Vol.103, pp. 280-292

Vol.132, No.7, pp.071014

*of Biomechanical Engineering*, Vol.102, pp.73-84

Behavior of Normal, Proteoglycan Depleted and Collagen Degraded Articular

Behavior of Articular Cartilage and Its Associated Electromechanical Properties.

Composition on Dynamic Visco-Elastic Property of Cartilaginous Tissue: Criteria for Classification between Hyaline Cartilage and Fibrocartilage Based on

Assessment of Fixed Charge Density in Regenerated Cartilage by Gd-DTPA - Enhanced MR Imaging. *Magnetic Resonance and Medical Science*, Vol.5, No.2, pp. 73-78

Noninvasive Evaluation of Biophysical Properties of Tissue-Engineered Cartilage

Fixed-charge Density in Tissue-Engineered Cartilage by Quantitative MRI and Relationship with Biomechanical Properties. *Jounal of Biomechanical Engineering*,

Relaxation of Articular Cartilage in Compression? Theory and Experiments. *Journal* 

Helminen, H.J. & Jurvelin, J.S. (2001). T2 Relaxation Reveals Spatial Collagen Architecture in Articular Cartilage: a Comparative Quantitative MRI and Polarized

Jurvelin, J.S. (2000). Quantitative MR Microscopy of Enzymatically Degraded

J.S. (2007). Estimation of Mechanical Properties of Articular Cartilage with MRI dGEMRIC, T2 and T1 Imaging in Different Species with Variable Stages of

R.G. (1998). Cartilage Formation in a Hollow Fiber Bioreactor Studied by Proton

Tissue to Biochemical Agents as Studied by Proton Magnetic Resonance

Noninvasive Assessment of Glycosaminoglycan Production in Injectable Tissue-Engineered Cartilage Constructs Using Magnetic Resonance Imaging. *Tissue* 

Distribution Patterns Inside Bovine Articular Cartilage as Visualized by 1H

Articular Cartilage: A Microscopic MRI (microMRI) Study. *Magnetic Resonance in* 

Cutaneous wounds, ulcers or burns, as a result of external damage, such as intensive solar exposure or accidents, occur in high numbers and may cause ever lasting traumas. In some cases, the wounds are very painful and impair an individual's daily life in terms of health, and social interaction with his or her surroundings, but also in terms of mobility and ability to work or to perform leisure activities. Furthermore, the aesthetic feeling of patients is highly improved by the treatment of burns and other wounds since otherwise the wounds can make a detrimental impact on the appearance of the patient's skin. Besides aesthetic appearance, also the functioning of the damaged skin is different since infections are possibly more likely to occur near the (scar tissue of the) burn. Furthermore, the conception of temperature and pain by a patient is also altered by the burn, ulcer or cutaneous wound. Therefore, the treatment of these defects, which occurs in very high numbers, is important and hence it is crucial to improve treatments and to make treatments more efficient.

Nowadays many different treatments, such as, in more serious cases, skin transplantation, injection of cell cultures and the implant of artificial skin, are commonly carried out by plastical surgeons in order to relieve the patient and to minimize the aesthetic effects of scar tissue as much as possible. In order to be able to improve medical treatment of scars, it is important to understand the biological mechanism behind scar formation in terms of plastic deformation and a regenerated excess of collageneous matrix. More information about the treatment of scars and the biology behind the scar formation can, among others, be found in Bloemen (2011); van der Veer (2009). The deformations are caused by the contractile mechanism due to the pulling behavior of the (myo-)fibroblasts. Dallon (2008) experimentally found that the nature of contraction is predominantly determined by the following two mechanisms:


These mechanisms were simulated by Dallon (2010) by the use of a spring model. Since Dallon's two-dimensional model tracks the forces exerted by each individual philapod or

Contraction and Angiogenesis 3

A Mathematical Model for Wound Contraction and Angiogenesis 491

Fig. 1. A schematic of the events during wound healing. The dermis and epidermis are

in this process on contraction and healing, as well as values such that healing of burns is optimized with respect to a minimal level of wound contraction and maximum healing speed. These results are helpful to surgeons to design innovative minimal invasive treatments that are applicable to patients. Furthermore, the mathematical analysis is helpful to make the

In this manuscript, we address two partial processes during the proliferative phase: wound contraction and angiogenesis. When the tissue is provided with enough oxygen and nutrients the process of wound closure starts. Cells in the epidermis, which mainly consist of keratinocytes, start regenerating the upperlayer of the wound. Usually the skin can not be replaced fully and some scars are left where the wound was located, due to excess of regenerated collagen. The second and third stage of wound healing do not take place at the same location in the wound. The former is located in the dermis, the latter is limited to the epidermis. The epidermis and the dermis consist of different type of cells and are separated

In this manuscript, we first present a selection of the currently available mathematical models that seek to describe the biological processes of wound healing as well as possible. The healing process is very complex and many factors contribute to it, therefore simplifications have to be made. The main topics in this thesis are combining two models for angiogenesis and coupling models for the different stages of wound healing. Currently models exist for angiogenesis and dermal regeneration (fibroblasts and collagen) separately, however, only scarcely have there been attempts to couple the models. This is vital as the various stages of wound healing overlap and hence influence each other. The models that we consider in this study are all formulated in terms of partial differential equations (PDEs) that are based on conservation principles. Such coupled models could give more insights into how the process of wound

illustrated. The picture was taken with permission from http://www.bioscience.org/2006/v11/af/1843/figures.htm

by a so-called basal membrane, see also Figure 1.

treatment patient specific.

lamellipodium of the individual (myo)fibroblast, the formalism becomes very expensive from a computational point of view when applied to burns of realistic dimensions with millions of fibroblasts and collageneous fibres. However, his observations and calculations are very useful for a thorough understanding of the implications of the fundamental processes and for the calibration of continuum models that are based on partial differential equations.

The severity of a wound is determined by the thickness of the damaged skin layer of the patient. A minor injury concerns the damage on the corneum, which is replaced with the patient's tissues relatively quickly. This type of injury is referred to as a first degree skin burn or wound. A more severe damage concerns the impairment of the epidermis. In this case reepithialialization has to repair the epidermis by migration and proliferation of keratinocytes. This damage is classified as a second degree burn or wound, but this type of wound will heal without many problems for healthy patients. Contraction does not take place and hence the burn or wound will recover entirely without any plastic deformations or significant scars. Third degree burns or wounds pose a more serious damage in the sense that also (part of) the dermis, and even the subcutis could be disrupted. Here the fibroblast rich dermis, consisting of the collageneous fibrous tissue has to be repaired. This takes place by a sequence of signaling processes to initiate proliferation and movement of fibroblasts and the restoration of the vascular network. As the fibroblasts produce an excess of collageneous fibres, on which they, and the myofibroblasts to a larger extent, exert contractile forces, wound contraction takes place. The process of wound contraction is a useful mechanism for a rapid minimization of exposure of the underlying tissues to hazardous external environments, when the wound is caused by mechanical damage. In particular, in skins of rabbits, percentages of up to 80 percent of the initial wound area have been reported in various experimental studies. In humans, experimental studies evidence that wound contraction is much less significant, being in the order of 5–10 percent. However, in the case of healing of burns, this contraction phenomenon is undesirable, as it gives rise to significant deformations, which can be plastic and furthermore, the patient is left with an excess of collageneous fibres and unpleasantly looking scars. Remodeling is a very slow process which cannot remove all wrinkles as a result of the plastic deformations due to wound contraction.

In order to be able to improve surgical treatments, in terms of effectiveness and minimization of invasion into the patient, it is crucial to know which behavior of the (myo-)fibroblasts cause the plastic deformation of the skin. In the case of a third degree burn or wound, fibroblasts enter the wound area during the proliferative stages of the healing process of the dermis. Subsequently, the fibroblasts start to proliferate and to produce extracellular matrix, referred to as fibrous tissue. Furthermore, as a result of exposure to high strains, fibroblasts differentiate into myofibroblasts, which are considered as weak muscle cells. Myofibroblasts increase the measure of contraction of the wound area of the dermis. This differentiation process also takes place under the influence of a growth factor TF-*β*, which is produced by the fibroblasts and myofibroblasts depending on the amount of extracellular matrix present. The growth factor also stimulates the production of fibroblasts and its regeneration of collagen. Next to the enhancement of the production of several biological entities, the growth factor increases the chemotactic transport of the fibroblasts towards the wound area. To get more quantitative insight into the influences of the growth factors and other agents on the amount of wound contraction caused by the pulling forces of the (myo)fibroblasts, detailed mathematical models with parameter sensitivity analysis are indispensable. This sensitivity analysis can give insight into the quantification of the influence of all biological parameters involved 2 Will-be-set-by-IN-TECH

lamellipodium of the individual (myo)fibroblast, the formalism becomes very expensive from a computational point of view when applied to burns of realistic dimensions with millions of fibroblasts and collageneous fibres. However, his observations and calculations are very useful for a thorough understanding of the implications of the fundamental processes and for

The severity of a wound is determined by the thickness of the damaged skin layer of the patient. A minor injury concerns the damage on the corneum, which is replaced with the patient's tissues relatively quickly. This type of injury is referred to as a first degree skin burn or wound. A more severe damage concerns the impairment of the epidermis. In this case reepithialialization has to repair the epidermis by migration and proliferation of keratinocytes. This damage is classified as a second degree burn or wound, but this type of wound will heal without many problems for healthy patients. Contraction does not take place and hence the burn or wound will recover entirely without any plastic deformations or significant scars. Third degree burns or wounds pose a more serious damage in the sense that also (part of) the dermis, and even the subcutis could be disrupted. Here the fibroblast rich dermis, consisting of the collageneous fibrous tissue has to be repaired. This takes place by a sequence of signaling processes to initiate proliferation and movement of fibroblasts and the restoration of the vascular network. As the fibroblasts produce an excess of collageneous fibres, on which they, and the myofibroblasts to a larger extent, exert contractile forces, wound contraction takes place. The process of wound contraction is a useful mechanism for a rapid minimization of exposure of the underlying tissues to hazardous external environments, when the wound is caused by mechanical damage. In particular, in skins of rabbits, percentages of up to 80 percent of the initial wound area have been reported in various experimental studies. In humans, experimental studies evidence that wound contraction is much less significant, being in the order of 5–10 percent. However, in the case of healing of burns, this contraction phenomenon is undesirable, as it gives rise to significant deformations, which can be plastic and furthermore, the patient is left with an excess of collageneous fibres and unpleasantly looking scars. Remodeling is a very slow process which cannot remove all wrinkles as a result

In order to be able to improve surgical treatments, in terms of effectiveness and minimization of invasion into the patient, it is crucial to know which behavior of the (myo-)fibroblasts cause the plastic deformation of the skin. In the case of a third degree burn or wound, fibroblasts enter the wound area during the proliferative stages of the healing process of the dermis. Subsequently, the fibroblasts start to proliferate and to produce extracellular matrix, referred to as fibrous tissue. Furthermore, as a result of exposure to high strains, fibroblasts differentiate into myofibroblasts, which are considered as weak muscle cells. Myofibroblasts increase the measure of contraction of the wound area of the dermis. This differentiation process also takes place under the influence of a growth factor TF-*β*, which is produced by the fibroblasts and myofibroblasts depending on the amount of extracellular matrix present. The growth factor also stimulates the production of fibroblasts and its regeneration of collagen. Next to the enhancement of the production of several biological entities, the growth factor increases the chemotactic transport of the fibroblasts towards the wound area. To get more quantitative insight into the influences of the growth factors and other agents on the amount of wound contraction caused by the pulling forces of the (myo)fibroblasts, detailed mathematical models with parameter sensitivity analysis are indispensable. This sensitivity analysis can give insight into the quantification of the influence of all biological parameters involved

the calibration of continuum models that are based on partial differential equations.

of the plastic deformations due to wound contraction.

Fig. 1. A schematic of the events during wound healing. The dermis and epidermis are illustrated. The picture was taken with permission from http://www.bioscience.org/2006/v11/af/1843/figures.htm

in this process on contraction and healing, as well as values such that healing of burns is optimized with respect to a minimal level of wound contraction and maximum healing speed. These results are helpful to surgeons to design innovative minimal invasive treatments that are applicable to patients. Furthermore, the mathematical analysis is helpful to make the treatment patient specific.

In this manuscript, we address two partial processes during the proliferative phase: wound contraction and angiogenesis. When the tissue is provided with enough oxygen and nutrients the process of wound closure starts. Cells in the epidermis, which mainly consist of keratinocytes, start regenerating the upperlayer of the wound. Usually the skin can not be replaced fully and some scars are left where the wound was located, due to excess of regenerated collagen. The second and third stage of wound healing do not take place at the same location in the wound. The former is located in the dermis, the latter is limited to the epidermis. The epidermis and the dermis consist of different type of cells and are separated by a so-called basal membrane, see also Figure 1.

In this manuscript, we first present a selection of the currently available mathematical models that seek to describe the biological processes of wound healing as well as possible. The healing process is very complex and many factors contribute to it, therefore simplifications have to be made. The main topics in this thesis are combining two models for angiogenesis and coupling models for the different stages of wound healing. Currently models exist for angiogenesis and dermal regeneration (fibroblasts and collagen) separately, however, only scarcely have there been attempts to couple the models. This is vital as the various stages of wound healing overlap and hence influence each other. The models that we consider in this study are all formulated in terms of partial differential equations (PDEs) that are based on conservation principles. Such coupled models could give more insights into how the process of wound

Contraction and Angiogenesis 5

A Mathematical Model for Wound Contraction and Angiogenesis 493

During the wound contraction stage fibroblasts (connective tissue cells) invade the wound site and contract the extracellular matrix (ECM). The contraction decreases the area of contact between the wound and its surroundings, thereby reducing the chance of contamination and infection. Furthermore this process is vital in assuring that new blood vessels can be formed in the wound during angiogenesis, since the fibroblasts invading the wound form the tissue in which the new capillaries can grow. The wound contraction stage is limited to the dermis,

We use the model for contraction due to Javierre (2009), which deals with the presence of myofibroblasts. Myofibroblasts are a kind of weak muscle cells. They are nonmotile cells that differentiate from fibroblasts and transmit and amplify the traction forces generated by the fibroblasts, Vermolen (2009). Secondly, the model incorporates the effects of a growth factor

*u*fib − *D*fib∇*u*fib +

where *c*ecm and *u*myo respectively denote the growth factor and myofibroblast concentration. The first term in the left-hand side corresponds to the total accumulation, the second term follows from passive convection due to deformation of the tissue. The third and fourth terms of the left-hand side account for random walk and chemotaxis (movement towards the gradient of the growth factor). In the right-hand side of the above equation, we respectively have the logistic proliferation term up to an equilibrium, in which proliferation is enhanced by the presence of the growth factor, the differentiation term to myofibroblasts under presence of the growth factor, back differiation from myofibroblasts and finally a term dealing with cell death. The cell death rate is denoted by *d*fib, the myofibroblast to fibroblast differentiation rate

*Ck* are known constants that monitor the growth factor's influence on the contraction process

The production term, first on the right hand side, now also incorporates growth factor stimulated proliferation. The other three terms on the right hand side respectively account

The PDE for the myofibroblast concentration *u*myo is similar to equation (1). However, since myofibroblast are nonmotile cells, they will only move due to passive convection. The

*λ*fib +

+

*λ*0 fib*c*ecm *<sup>C</sup>*1/2 <sup>+</sup> *<sup>c</sup>*ecm

*k*1*c*ecm *Ck* + *c*ecm *<sup>u</sup>*myo

<sup>1</sup> <sup>−</sup> *<sup>u</sup>*myo *K* 

*u*fib − *k*2*u*myo − *d*myo*u*myo, (2)

<sup>−</sup> *<sup>k</sup>*1*c*ecm *Ck* + *c*ecm

<sup>1</sup> <sup>−</sup> *<sup>u</sup>*fib *K* 

by *k*<sup>2</sup> and the fibroblast to myofibroblast differentiation rate by *k*1. Furthermore *λ*<sup>0</sup>

<sup>=</sup> *<sup>ε</sup>*myo

and *K* is a parameter that regulates the equilibrium concentration.

for differentiation to and from myofibroblasts and cell death.

*<sup>u</sup>*myo

myofibroblast concentration thus obeys

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · *<sup>∂</sup>***<sup>u</sup>**

*∂t*

*∂u*myo

*a*fib

(*b*fib <sup>+</sup> *<sup>c</sup>*ecm)<sup>2</sup> *<sup>u</sup>*fib∇*c*ecm

=

fib, *C*1/2 and

*u*fib + *k*2*u*myo − *d*fib*u*fib, (1)

however, the contraction of the ECM also effects the tissue in the epidermis.

The equation concerning the fibroblast concentration *u*fib becomes

*∂t*

*<sup>u</sup>*fib

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · *<sup>∂</sup>***<sup>u</sup>**

**2.1 The wound contraction models**

that triggers wound contraction.

*λ*fib +

*∂u*fib

*λ*0 fib*c*ecm *<sup>C</sup>*1/2 <sup>+</sup> *<sup>c</sup>*ecm

healing evolves. These insights might lead to treatments that reduce healing time, e.g. the use of certain hormones to speed up the healing process. Also scars and other deformations due to incomplete healing might be prevented or reduced. This issue is especially crucial for the treatment of burns, where hypotrophic scars may result after injury.

A lot of simulations have been done in literature to give more insight on how these models behave. Also the dependence of the models on certain parameters is investigated in the present paper. Finally recommendations are made for future research in the topic of mathematically describing wound healing. Many studies are based on PDEs, for instance the work by Murray (2004); Olsen (1995); Maggelakis (2003); Gaffney (2002); Sherratt (1991); Adam (1999); Javierre (2009); Vermolen (2009; 2010; 2011); Schugart (2008); Xue (2009). Many other modeling studies are based on cell-based or Monte-Carlo methods, such as the cellular Potts models used by Glazier (1992); Merks (2009); Plank (2004) to mention a few. This class of model is lattice-based and minimizes a virtual energy functional. Recently, a semi-stochastic continuous cell-based model was formulated by Vermolen-Gefen (2011).

The present manuscript falls within the class of papers, in which one attempts to combine existing mathematical models for various subprocesses occuring during wound healing. Some recent work by this team are Vermolen (2010), and Vermolen (2011), in which the most simplified models for wound contraction, closure and angiogenesis are coupled. The present model is based on a combination of some more advanced models for wound contraction and angiogenesis and gives a more sophisticated formalism for dermal regeneration. The paper is organized as follows. In Section 2, we present several models for processes of angiogenesis and wound contraction. Here, we also present some results of simulations. This is an innovation with respect to other studies. In Section 3, we consider the coupling of the models, including simulations. In Section 4, we describe the numerical methods. Finally, some conclusions are drawn in Section 5.

The character of the current manuscript is rather informal and descriptive about the mathematical concepts used in the present study. The level of abstract mathematics has been reduced tremendously. Further, we note that the simulations shown in this manuscript are still in a preliminary state.

#### **2. Current models**

The mathematical model for the proliferative stage of wound healing is usually separated in three distinct parts representing three stages of wound healing. These three stages are wound contraction, angiogenesis and wound closure. Note that the inflammation stage mentioned in Section 1 is not taken into account. This is due to the fact that inflammation only contains the damage and only after the inflammation stage is finished the real healing process starts.

In this chapter we present some of the currently available models on the above mentioned wound healing stages. In Section 2.1, we first present the model on wound contraction. Next in Section 2.2, a model for angiogenesis is presented. This model is a combination of two accepted models that describe a specific feature in the angiogenesis process. The two models take a very different approach on how to model the growth of new blood vessels in the wound. In Section 2.3, a study that attempts to combine models of dermal regeneration and angiogenesis is briefly discussed. In all the models that we deal with, we consider a bounded simply connected domain <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>**2. The boundary is denoted by *<sup>∂</sup>*Ω.

#### **2.1 The wound contraction models**

4 Will-be-set-by-IN-TECH

healing evolves. These insights might lead to treatments that reduce healing time, e.g. the use of certain hormones to speed up the healing process. Also scars and other deformations due to incomplete healing might be prevented or reduced. This issue is especially crucial for the

A lot of simulations have been done in literature to give more insight on how these models behave. Also the dependence of the models on certain parameters is investigated in the present paper. Finally recommendations are made for future research in the topic of mathematically describing wound healing. Many studies are based on PDEs, for instance the work by Murray (2004); Olsen (1995); Maggelakis (2003); Gaffney (2002); Sherratt (1991); Adam (1999); Javierre (2009); Vermolen (2009; 2010; 2011); Schugart (2008); Xue (2009). Many other modeling studies are based on cell-based or Monte-Carlo methods, such as the cellular Potts models used by Glazier (1992); Merks (2009); Plank (2004) to mention a few. This class of model is lattice-based and minimizes a virtual energy functional. Recently, a semi-stochastic

The present manuscript falls within the class of papers, in which one attempts to combine existing mathematical models for various subprocesses occuring during wound healing. Some recent work by this team are Vermolen (2010), and Vermolen (2011), in which the most simplified models for wound contraction, closure and angiogenesis are coupled. The present model is based on a combination of some more advanced models for wound contraction and angiogenesis and gives a more sophisticated formalism for dermal regeneration. The paper is organized as follows. In Section 2, we present several models for processes of angiogenesis and wound contraction. Here, we also present some results of simulations. This is an innovation with respect to other studies. In Section 3, we consider the coupling of the models, including simulations. In Section 4, we describe the numerical methods. Finally, some

The character of the current manuscript is rather informal and descriptive about the mathematical concepts used in the present study. The level of abstract mathematics has been reduced tremendously. Further, we note that the simulations shown in this manuscript are

The mathematical model for the proliferative stage of wound healing is usually separated in three distinct parts representing three stages of wound healing. These three stages are wound contraction, angiogenesis and wound closure. Note that the inflammation stage mentioned in Section 1 is not taken into account. This is due to the fact that inflammation only contains the damage and only after the inflammation stage is finished the real healing process starts.

In this chapter we present some of the currently available models on the above mentioned wound healing stages. In Section 2.1, we first present the model on wound contraction. Next in Section 2.2, a model for angiogenesis is presented. This model is a combination of two accepted models that describe a specific feature in the angiogenesis process. The two models take a very different approach on how to model the growth of new blood vessels in the wound. In Section 2.3, a study that attempts to combine models of dermal regeneration and angiogenesis is briefly discussed. In all the models that we deal with, we consider a

bounded simply connected domain <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>**2. The boundary is denoted by *<sup>∂</sup>*Ω.

treatment of burns, where hypotrophic scars may result after injury.

continuous cell-based model was formulated by Vermolen-Gefen (2011).

conclusions are drawn in Section 5.

still in a preliminary state.

**2. Current models**

During the wound contraction stage fibroblasts (connective tissue cells) invade the wound site and contract the extracellular matrix (ECM). The contraction decreases the area of contact between the wound and its surroundings, thereby reducing the chance of contamination and infection. Furthermore this process is vital in assuring that new blood vessels can be formed in the wound during angiogenesis, since the fibroblasts invading the wound form the tissue in which the new capillaries can grow. The wound contraction stage is limited to the dermis, however, the contraction of the ECM also effects the tissue in the epidermis.

We use the model for contraction due to Javierre (2009), which deals with the presence of myofibroblasts. Myofibroblasts are a kind of weak muscle cells. They are nonmotile cells that differentiate from fibroblasts and transmit and amplify the traction forces generated by the fibroblasts, Vermolen (2009). Secondly, the model incorporates the effects of a growth factor that triggers wound contraction.

The equation concerning the fibroblast concentration *u*fib becomes

$$\frac{\partial u\_{\rm fib}}{\partial t} + \nabla \cdot \left( \frac{\partial \mathbf{u}}{\partial t} u\_{\rm fib} - D\_{\rm fib} \nabla u\_{\rm fib} + \frac{a\_{\rm fib}}{(b\_{\rm fib} + c\_{\rm ecm})^2} u\_{\rm fib} \nabla c\_{\rm ecm} \right) =$$

$$\left( \lambda\_{\rm fib} + \frac{\lambda\_{\rm fib}^0 c\_{\rm ecm}}{\mathbb{C}\_{1/2} + c\_{\rm ecm}} \right) u\_{\rm fib} \left( 1 - \frac{u\_{\rm fib}}{K} \right) - \frac{k\_1 c\_{\rm ecm}}{\mathbb{C}\_k + c\_{\rm ecm}} u\_{\rm fib} + k\_2 u\_{\rm proj} - d\_{\rm fib} u\_{\rm fib} \tag{1}$$

where *c*ecm and *u*myo respectively denote the growth factor and myofibroblast concentration. The first term in the left-hand side corresponds to the total accumulation, the second term follows from passive convection due to deformation of the tissue. The third and fourth terms of the left-hand side account for random walk and chemotaxis (movement towards the gradient of the growth factor). In the right-hand side of the above equation, we respectively have the logistic proliferation term up to an equilibrium, in which proliferation is enhanced by the presence of the growth factor, the differentiation term to myofibroblasts under presence of the growth factor, back differiation from myofibroblasts and finally a term dealing with cell death. The cell death rate is denoted by *d*fib, the myofibroblast to fibroblast differentiation rate by *k*<sup>2</sup> and the fibroblast to myofibroblast differentiation rate by *k*1. Furthermore *λ*<sup>0</sup> fib, *C*1/2 and *Ck* are known constants that monitor the growth factor's influence on the contraction process and *K* is a parameter that regulates the equilibrium concentration.

The production term, first on the right hand side, now also incorporates growth factor stimulated proliferation. The other three terms on the right hand side respectively account for differentiation to and from myofibroblasts and cell death.

The PDE for the myofibroblast concentration *u*myo is similar to equation (1). However, since myofibroblast are nonmotile cells, they will only move due to passive convection. The myofibroblast concentration thus obeys

$$\begin{aligned} \frac{\partial u\_{\rm myo}}{\partial t} + \nabla \cdot \left( \frac{\partial \mathbf{u}}{\partial t} u\_{\rm myo} \right) &= \varepsilon\_{\rm myo} \left( \lambda\_{\rm fib} + \frac{\lambda\_{\rm fib}^{0} \varepsilon\_{\rm ccm}}{\mathbb{C}\_{1/2} + \mathfrak{c}\_{\rm ccm}} \right) u\_{\rm myo} \left( 1 - \frac{u\_{\rm myo}}{K} \right) \\ &+ \frac{k\_{1} \varepsilon\_{\rm ccm}}{\mathbb{C}\_{k} + \mathfrak{c}\_{\rm ccm}} u\_{\rm fib} - k\_{2} u\_{\rm myo} - d\_{\rm myo} u\_{\rm myo} \,, \end{aligned} \tag{2}$$

Contraction and Angiogenesis 7

A Mathematical Model for Wound Contraction and Angiogenesis 495

The mechanical part of the wound contraction model is based on the linear viscoelastic

Here *σ* = *σ*ecm + *σ*cell, the stresstensor, accounts for the ECM related stress, *σ*ecm, and the cell

Here the first two terms on the right hand side respresent the viscous effects and the last term the elastic effects. If we let **u** = **u**(**x**, *t*) denote the displacement of the ECM, then the strain

Furthermore, in relation (6), **I** denotes the identity tensor and *μ*1, *μ*2, *E* and *ν* respectively

At time *t* = 0 it is assumed that there is no displacement of the ECM, i.e. **u**(**x**, 0) = **0**. Also we assume that **u** vanishes at the boundary far away from the wound, i.e. **u**(**x**, *t*) = **0** for **x** ∈ *∂*Ω. This can be justified by taking the computational domain Ω sufficiently large, so that

Since the myofibroblasts transmit and amplify the traction forces generated by the fibroblasts this is also visible in the cell traction term *σ*cell. For the model due to Olsen et al. this term is

*R*2

where *ξ* is a constant of proportionality and *R<sup>τ</sup>* quantifies how the cell traction depends on

In Javierre (2009) an extension of the model due to Olsen et al. is presented. Javierre et al. propose the mechanical stress to act as a factor that effects the differentiation from fibroblasts to myofibroblasts. They introduce an estimation of the mechanical stimulus that depends on

(*θ*<sup>1</sup> <sup>−</sup> *<sup>θ</sup>*) *<sup>χ</sup>*[*θ*1,*θ*∗](*θ*) + *<sup>K</sup>*act*p*max

1 + *ξu*myo

*K*act*θ*<sup>2</sup> − *p*max

*ρ*

The external forces acting on the tissue, **f**ext, are modelled similarly in all three models, i.e.

represent the dynamic and kinematic viscosity, Young's modulus and Poisson's ratio.

Here *ρ* = *ρ*(**x**, *t*) denotes the ECM density and *s* is the tethering elasticity coefficient.

*<sup>σ</sup>*cell <sup>=</sup> *<sup>τ</sup>u*fib

∇**u** + (∇**u**)

*E* 1 + *ν*  *�* +

> *T*

stress, *σ*cell. Furthermore **f**ext represents the external forces acting on the tissue. In all the below discussed models the ECM related stress tensor *σ*ecm is given as

> *∂θ ∂t* **I** +

*σ*ecm = *μ*<sup>1</sup>

the boundary effects can be ignored.

*∂� <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup>

tensor *�* and the dilation *θ* in equation (6) are respectively given by

*�* <sup>=</sup> <sup>1</sup> 2 

−∇· *σ* = **f**ext. (5)

*θ***I** 

*θ* = ∇ · **u**. (8)

**f**ext = −*sρ***u**. (9)

*<sup>τ</sup>* <sup>+</sup> *<sup>ρ</sup>*<sup>2</sup> **<sup>I</sup>**, (10)

(*θ*<sup>2</sup> − *θ*) *χ*[*θ*2,*θ*∗](*θ*) + *K*pas*θ*. (11)

. (6)

(7)

*ν* 1 − 2*ν*

equations, i.e.

and

given by

the ECM density.

the dilation *θ* = ∇ · **u** as *<sup>p</sup>*cell(*θ*) = *<sup>K</sup>*act*p*max

*K*act*θ*<sup>1</sup> − *p*max

where *ε*ecm is a proportionality constant, further *d*myo and *u*<sup>0</sup> myo denote the myofibroblasts death rate and the myofibroblast equilibrium concentration respectively. The terms on the left-hand side account for accumulation and passive convection due to deformation of the tissue. The right-hand side contains proliferation of myofibroblasts under presence of the growth factor, differentiation of fibroblasts to myofibroblasts, back differentiation to fibroblasts, and (programmed) cell death (apoptosis).

Both fibroblasts and myofibroblasts contribute to the production of the ECM, furthermore the production is chemically enhanced by the growth factor, Vermolen (2009). The PDE for the ECM density *ρ* then is given by

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\frac{\partial \mathbf{u}}{\partial t} \rho\right) = \left(\lambda\_{\rho} + \frac{\lambda\_{\rho}^{0} c\_{\text{ecm}}}{\mathbb{C}\_{\rho} + c\_{\text{ecm}}}\right) \frac{u\_{\text{fib}} + \eta\_{b} u\_{\text{myo}}}{R\_{\rho}^{2} + \rho^{2}} - d\_{\rho} (u\_{\text{fib}} + \eta\_{d} u\_{\text{myo}}) \rho\_{\text{}} \tag{3}$$

where *λρ* and *d<sup>ρ</sup>* are the ECM production and death rate respectively. The terms on the left-hand side describe accumulation and passive convection, and the right-hand side describes growth factor enhanced production of collagen, as well as decay of collagen. Furthermore *λ*<sup>0</sup> *<sup>ρ</sup>* and *C<sup>ρ</sup>* are known constants that monitor the growth factor's influence on the contraction process. Further, *R<sup>ρ</sup>* is a parameter that quantifies how the ECM production rate depends on the ECM density itself and *η<sup>b</sup>* and *η<sup>d</sup>* are proportionality constants.

The dynamics of the growth factor concentration *c*ecm are mainly determined by the fibroblasts and myofibroblasts as they produce the growth factor. Also the growth factor is motile, so it is subject to active convection. This leads to the following PDE for the growth factor concentration

$$\frac{\partial \mathbf{c}\_{\rm ecm}}{\partial t} + \nabla \cdot \left( \frac{\partial \mathbf{u}}{\partial t} \mathbf{c}\_{\rm ecm} - D\_{\rm c} \nabla c\_{\rm ecm} \right) = \frac{k\_{\rm c} (u\_{\rm fib} + \zeta u\_{\rm myo}) c\_{\rm ecm}}{\Gamma + c\_{\rm ecm}} - d\_{\rm c} c\_{\rm ecm}.\tag{4}$$

Here the second and third term on the left hand side account for passive convection and ordinary diffusion. The terms on the right hand side account for growth factor production and growth factor decay respectively. Furthermore *Dc* is the growth factor diffusion coefficient, *kc* denotes the growth factor production rate and *dc* the natural decay rate. Also Γ is a parameter that quantifies how the growth factor production rate depends on the growth factor concentration itself and *ζ* is a proportionality constant.

The initial and boundary conditions for the fibroblast concentration and the ECM density remain the same as in the model due to Tranquillo. For the myofibroblast and growth factor concentration the following initial conditions are imposed

$$
\mu\_{\rm myo}(\mathbf{x},0) = 0, c\_{\rm ecm}(\mathbf{x},0) = c\_{\rm ecm}^0
$$

for **x** ∈ Ω*<sup>w</sup>* and

$$u\_{\rm myo}(\mathbf{x},0) = 0, c\_{\rm ecm}(\mathbf{x},0) = 0$$

for **<sup>x</sup>** <sup>∈</sup> <sup>Ω</sup>*u*. Here *<sup>c</sup>*<sup>0</sup> ecm denotes the growth factor equilibrium concentration. Furthermore the growth factor concentration satisfy a no flux boundary condition on all boundaries. For the myofibroblast concentration we can apply the same reasoning as for the ECM density *ρ* and hence no boundary conditions have to be imposed.

The mechanical part of the wound contraction model is based on the linear viscoelastic equations, i.e.

$$-\nabla \cdot \boldsymbol{\sigma} = \mathbf{f}\_{\text{ext}}.\tag{5}$$

Here *σ* = *σ*ecm + *σ*cell, the stresstensor, accounts for the ECM related stress, *σ*ecm, and the cell stress, *σ*cell. Furthermore **f**ext represents the external forces acting on the tissue.

In all the below discussed models the ECM related stress tensor *σ*ecm is given as

$$
\sigma\_{\rm ecm} = \mu\_1 \frac{\partial \epsilon}{\partial t} + \mu\_2 \frac{\partial \theta}{\partial t} \mathbf{I} + \frac{E}{1+\nu} \left( \varepsilon + \frac{\nu}{1-2\nu} \theta \mathbf{I} \right). \tag{6}
$$

Here the first two terms on the right hand side respresent the viscous effects and the last term the elastic effects. If we let **u** = **u**(**x**, *t*) denote the displacement of the ECM, then the strain tensor *�* and the dilation *θ* in equation (6) are respectively given by

$$\boldsymbol{\epsilon} = \frac{1}{2} \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right) \tag{7}$$

and

6 Will-be-set-by-IN-TECH

death rate and the myofibroblast equilibrium concentration respectively. The terms on the left-hand side account for accumulation and passive convection due to deformation of the tissue. The right-hand side contains proliferation of myofibroblasts under presence of the growth factor, differentiation of fibroblasts to myofibroblasts, back differentiation to

Both fibroblasts and myofibroblasts contribute to the production of the ECM, furthermore the production is chemically enhanced by the growth factor, Vermolen (2009). The PDE for the

where *λρ* and *d<sup>ρ</sup>* are the ECM production and death rate respectively. The terms on the left-hand side describe accumulation and passive convection, and the right-hand side describes growth factor enhanced production of collagen, as well as decay of collagen.

the contraction process. Further, *R<sup>ρ</sup>* is a parameter that quantifies how the ECM production

The dynamics of the growth factor concentration *c*ecm are mainly determined by the fibroblasts and myofibroblasts as they produce the growth factor. Also the growth factor is motile, so it is subject to active convection. This leads to the following PDE for the growth factor

Here the second and third term on the left hand side account for passive convection and ordinary diffusion. The terms on the right hand side account for growth factor production and growth factor decay respectively. Furthermore *Dc* is the growth factor diffusion coefficient, *kc* denotes the growth factor production rate and *dc* the natural decay rate. Also Γ is a parameter that quantifies how the growth factor production rate depends on the growth factor

The initial and boundary conditions for the fibroblast concentration and the ECM density remain the same as in the model due to Tranquillo. For the myofibroblast and growth factor

*u*myo(**x**, 0) = 0, *c*ecm(**x**, 0) = *c*<sup>0</sup>

*u*myo(**x**, 0) = 0, *c*ecm(**x**, 0) = 0

growth factor concentration satisfy a no flux boundary condition on all boundaries. For the myofibroblast concentration we can apply the same reasoning as for the ECM density *ρ* and

rate depends on the ECM density itself and *η<sup>b</sup>* and *η<sup>d</sup>* are proportionality constants.

*c*ecm − *Dc*∇*c*ecm

 *<sup>u</sup>*fib <sup>+</sup> *<sup>η</sup>bu*myo *R*2

*<sup>ρ</sup>* and *C<sup>ρ</sup>* are known constants that monitor the growth factor's influence on

<sup>=</sup> *kc* (*u*fib <sup>+</sup> *<sup>ζ</sup>u*myo)*c*ecm Γ + *c*ecm

ecm

ecm denotes the growth factor equilibrium concentration. Furthermore the

*λ*0 *<sup>ρ</sup>c*ecm *C<sup>ρ</sup>* + *c*ecm myo denote the myofibroblasts

*<sup>ρ</sup>* <sup>+</sup> *<sup>ρ</sup>*<sup>2</sup> <sup>−</sup> *<sup>d</sup>ρ*(*u*fib <sup>+</sup> *<sup>η</sup>du*myo)*ρ*, (3)

− *dcc*ecm. (4)

where *ε*ecm is a proportionality constant, further *d*myo and *u*<sup>0</sup>

fibroblasts, and (programmed) cell death (apoptosis).

*∂***u** *∂t*

concentration itself and *ζ* is a proportionality constant.

concentration the following initial conditions are imposed

hence no boundary conditions have to be imposed.

ECM density *ρ* then is given by

 *∂***u** *∂t ρ* = *λρ* +

*∂ρ <sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ ·

Furthermore *λ*<sup>0</sup>

concentration

for **x** ∈ Ω*<sup>w</sup>* and

for **<sup>x</sup>** <sup>∈</sup> <sup>Ω</sup>*u*. Here *<sup>c</sup>*<sup>0</sup>

*∂c*ecm *<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ ·

$$
\theta = \nabla \cdot \mathbf{u}.\tag{8}
$$

Furthermore, in relation (6), **I** denotes the identity tensor and *μ*1, *μ*2, *E* and *ν* respectively represent the dynamic and kinematic viscosity, Young's modulus and Poisson's ratio.

The external forces acting on the tissue, **f**ext, are modelled similarly in all three models, i.e.

$$\mathbf{f}\_{\text{ext}} = -s\rho \mathbf{u}.\tag{9}$$

Here *ρ* = *ρ*(**x**, *t*) denotes the ECM density and *s* is the tethering elasticity coefficient.

At time *t* = 0 it is assumed that there is no displacement of the ECM, i.e. **u**(**x**, 0) = **0**. Also we assume that **u** vanishes at the boundary far away from the wound, i.e. **u**(**x**, *t*) = **0** for **x** ∈ *∂*Ω. This can be justified by taking the computational domain Ω sufficiently large, so that the boundary effects can be ignored.

Since the myofibroblasts transmit and amplify the traction forces generated by the fibroblasts this is also visible in the cell traction term *σ*cell. For the model due to Olsen et al. this term is given by

$$
\sigma\_{\rm cell} = \frac{\tau \mu\_{\rm fib} \left( 1 + \xi \mu\_{\rm myo} \right) \rho}{R\_{\tau}^2 + \rho^2} \mathbf{I}\_{\prime} \tag{10}
$$

where *ξ* is a constant of proportionality and *R<sup>τ</sup>* quantifies how the cell traction depends on the ECM density.

In Javierre (2009) an extension of the model due to Olsen et al. is presented. Javierre et al. propose the mechanical stress to act as a factor that effects the differentiation from fibroblasts to myofibroblasts. They introduce an estimation of the mechanical stimulus that depends on the dilation *θ* = ∇ · **u** as

$$p\_{\rm cell}(\theta) = \frac{K\_{\rm act} p\_{\rm max}}{K\_{\rm act} \theta\_1 - p\_{\rm max}} \left(\theta\_1 - \theta\right) \chi\_{[\theta\_1 \theta^\*]}(\theta) + \frac{K\_{\rm act} p\_{\rm max}}{K\_{\rm act} \theta\_2 - p\_{\rm max}} \left(\theta\_2 - \theta\right) \chi\_{[\theta\_2 \theta^\*]}(\theta) + K\_{\rm pass} \theta. \tag{11}$$

Contraction and Angiogenesis 9

A Mathematical Model for Wound Contraction and Angiogenesis 497

−1 −0.8 −0.6 −0.4 −0.2 <sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> −0.5

Fig. 2. The mechanical stimulus *p*cell as a function of the dilation *θ*, computed with values

0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 3. Normalized fibroblast (left) and myofibroblast (right) densities two days after injury.

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

y

θ

0

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

x

All input data have been given in the appendix.

0.5

1

1.5

Pcell

from Javierre (2009).

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

y

2

2.5

<sup>3</sup> x 10−4

Here *χ* denotes the indicator function, i.e.

$$\chi\_I(\theta) = \begin{cases} 1, \text{ if } \theta \in I, \\ 0, \text{ else.} \end{cases}$$

and the first two terms on the right hand side account for the contractile stress generated internally by the myosin machinery and transmitted through the actin bundles, Javierre (2009). The third term on the right hand side establishes the contractile stress supported by the passive resistance of the cell. See also Figure 2 for a plot of *p*cell(*θ*).

Furthermore, in equation (11), the compression and traction strain limits are respectively denoted by *θ*<sup>1</sup> and *θ*2, *p*max respresents the maximal contractile force exerted by the actomyosin machinery and *K*max and *K*pas the volumetric stiffness moduli of the active and passive components of the cell. Also the parameter *θ*∗ can be computed from *K*act and *p*max as *θ*<sup>∗</sup> = *<sup>p</sup>*max *<sup>K</sup>*act , Javierre (2009).

To incorporate the effects of the mechanical stimulus on the fibroblast to myofibroblast differentiation an extra factor is found before the differentiation term (the second term on the right hand side) in equation (1). The fibroblast to myofibroblast differentiation term changes to

$$\frac{p\_{\text{cell}}(\theta)}{\pi\_d + p\_{\text{cell}}(\theta)} \frac{k\_1 c\_{\text{ecm}}}{\mathcal{C}\_k + c\_{\text{ecm}}} u\_{\text{fib}} \tag{12}$$

where *τ<sup>d</sup>* is a parameter that quantifies how the differentiation rate depends on the mechanical stimulus. Note that this term is also present in the PDE for the myofibroblast concentration and that thus the second term on the right hand side of equation (2) also changes to the above.

The mechanical stimulus also effects the the cell stresses and thus the cell traction term **œ**cell. In Javierre (2009) it is assumed that **œ**cell depends linearly on *p*cell(*θ*) and thus that

$$
\sigma\_{\rm cell} = p\_{\rm cell}(\theta) \frac{u\_{\rm filb} \left(1 + \xi u\_{\rm myo}\right) \rho}{R\_{\tau}^2 + \rho^2} \mathbf{I}. \tag{13}
$$

In order to give more insight in the process of wound contraction we did some simulations with the model due to Javierre (2009). We will describe the numerical techniques in Section 4. Further, the input data can be found in the appendix. As a computational domain we use the unit square, i.e. Ω = {**x** = (*x*, *y*)| 0 ≤ *x*, *y* ≤ 1}, and the initial wound is given by Ω*<sup>w</sup>* = **<sup>x</sup>**| |**x**| ≤ <sup>1</sup> 2 . The results are given in Figures 3 to 5, where we show the solution two days after injury. The computations have been done using parameter values taken from Javierre (2009), see also the appendix.

In Figure 3 we show the fibroblast and myofibroblast concentration two days after injury. The myofibroblast concentration is highly concentrated around the wound edge, whereas the fibroblast have invaded the wound. Figure 4 shows the ECM density and the growth factor concentration two days after injury. We see that the ECM density is slightly elevated at the wound edge, which is to be expected since the wound is healing there. Furthermore the growth factor concentration has spread throughout the computational domain, but is still concentrated inside the wound. In Figure 5, the displacement pattern due to contractile forces exerted by fibroblasts and myofibroblasts is plotted. If we compare Figure 5 with Figure 3 we see that the ECM displacement is largest at places of high myofibroblast concentration. Here 8 Will-be-set-by-IN-TECH

and the first two terms on the right hand side account for the contractile stress generated internally by the myosin machinery and transmitted through the actin bundles, Javierre (2009). The third term on the right hand side establishes the contractile stress supported by

Furthermore, in equation (11), the compression and traction strain limits are respectively denoted by *θ*<sup>1</sup> and *θ*2, *p*max respresents the maximal contractile force exerted by the actomyosin machinery and *K*max and *K*pas the volumetric stiffness moduli of the active and passive components of the cell. Also the parameter *θ*∗ can be computed from *K*act and *p*max as

To incorporate the effects of the mechanical stimulus on the fibroblast to myofibroblast differentiation an extra factor is found before the differentiation term (the second term on the right hand side) in equation (1). The fibroblast to myofibroblast differentiation term changes

where *τ<sup>d</sup>* is a parameter that quantifies how the differentiation rate depends on the mechanical stimulus. Note that this term is also present in the PDE for the myofibroblast concentration and that thus the second term on the right hand side of equation (2) also changes to the above. The mechanical stimulus also effects the the cell stresses and thus the cell traction term **œ**cell.

> *u*fib

In order to give more insight in the process of wound contraction we did some simulations with the model due to Javierre (2009). We will describe the numerical techniques in Section 4. Further, the input data can be found in the appendix. As a computational domain we use the unit square, i.e. Ω = {**x** = (*x*, *y*)| 0 ≤ *x*, *y* ≤ 1}, and the initial wound is given by

two days after injury. The computations have been done using parameter values taken from

In Figure 3 we show the fibroblast and myofibroblast concentration two days after injury. The myofibroblast concentration is highly concentrated around the wound edge, whereas the fibroblast have invaded the wound. Figure 4 shows the ECM density and the growth factor concentration two days after injury. We see that the ECM density is slightly elevated at the wound edge, which is to be expected since the wound is healing there. Furthermore the growth factor concentration has spread throughout the computational domain, but is still concentrated inside the wound. In Figure 5, the displacement pattern due to contractile forces exerted by fibroblasts and myofibroblasts is plotted. If we compare Figure 5 with Figure 3 we see that the ECM displacement is largest at places of high myofibroblast concentration. Here

*k*1*c*ecm *Ck* + *c*ecm

1 + *ξu*myo

. The results are given in Figures 3 to 5, where we show the solution

*R*2

 *ρ*

*u*fib, (12)

*<sup>τ</sup>* <sup>+</sup> *<sup>ρ</sup>*<sup>2</sup> **<sup>I</sup>**. (13)

*p*cell(*θ*) *τ<sup>d</sup>* + *p*cell(*θ*)

In Javierre (2009) it is assumed that **œ**cell depends linearly on *p*cell(*θ*) and thus that

*σ*cell = *p*cell(*θ*)

 1, if *<sup>θ</sup>* <sup>∈</sup> *<sup>I</sup>*, 0, else.

*χI*(*θ*) =

the passive resistance of the cell. See also Figure 2 for a plot of *p*cell(*θ*).

Here *χ* denotes the indicator function, i.e.

*θ*<sup>∗</sup> = *<sup>p</sup>*max

to

Ω*<sup>w</sup>* =

**<sup>x</sup>**| |**x**| ≤ <sup>1</sup> 2 

Javierre (2009), see also the appendix.

*<sup>K</sup>*act , Javierre (2009).

Fig. 2. The mechanical stimulus *p*cell as a function of the dilation *θ*, computed with values from Javierre (2009).

x

Fig. 3. Normalized fibroblast (left) and myofibroblast (right) densities two days after injury. All input data have been given in the appendix.

x

x

Contraction and Angiogenesis 11

A Mathematical Model for Wound Contraction and Angiogenesis 499

Fibroblasts Myofibroblasts ECM Growth Factor

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> <sup>7</sup> <sup>8</sup> <sup>9</sup> <sup>10</sup> <sup>0</sup>

Time in days

remark that all the simulations in this subsection were obtained without any coupling with

In this section, we present the model for angiogenesis. These two models each take a very different approach on how to model this process. Where the model due to Maggelakis (2003) focusses on the relation between a lack of oxygen and capillary growth, the model due to Gaffney (2002), attempts to model the migration of endothelial cells into the wound. Both models adress an important aspect of angiogenesis, but on the other hand both also miss an

In this section we present a novel angiogenesis model based on the models of Maggelakis and Gaffney et al. In this model we attempt to unite the two approaches and thus create a model that is suitable for dealing with both aspects of angiogenesis. This gives a new model for angiogenesis. By combining these two aspects in one model we hope to create a more accurate model of the process of angiogenesis. As a basis we take the model of Gaffney et al. and extend it with a negative feedback mechanism for oxygen. This assures that both the endothelial cell migration and the oxygen shortage aspects are covered in this novel angiogenesis model.

−∇· (*D*1∇*u*tip + *D*2∇*u*end) = *f*(*u*tip, *u*end),

−∇· *λ*5(*D*1∇*u*tip + *D*2∇*u*end) = *g*(*u*tip, *u*end)

(14)

Fig. 6. Normalized solutions furthest in the wound, at **x** = **0**, in time.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

The model due to Gaffney et al. is given by

*∂u*tip *∂t*

*∂u*end *∂t*

the model for angiogenesis.

**2.2 Angiogenesis**

important aspect.

Fig. 4. Normalized ECM density (left) and growth factor concentration (right) two days after injury.

Fig. 5. Displacement of the ECM two days after injury.

x

the effect of the myofibroblasts, as they transmit and amplify the traction forces generated by the fibroblasts, can clearly be seen.

The normalized solutions in time at **x** = **0**, furthest in the wound, are shown in Figure 6. This gives an indication of the development of the solutions in time. We see that the fibroblast concentration increases towards it equilibrium, whereas the myofibroblast concentration first rises and then falls again. This is due to the fact that, as the wound heals, the myofibroblasts either differentiate back to fibroblasts again or undergo apoptosis and hence eventually will disappear completely. Also the growth factor concentration eventually goes to zero as the wound is healed. The ECM density slowly rises in time and will eventually reach its equilibrium, allthough it can take some time before the ECM is completely restored. We finally

Fig. 6. Normalized solutions furthest in the wound, at **x** = **0**, in time.

remark that all the simulations in this subsection were obtained without any coupling with the model for angiogenesis.

#### **2.2 Angiogenesis**

10 Will-be-set-by-IN-TECH

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 4. Normalized ECM density (left) and growth factor concentration (right) two days after

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

the effect of the myofibroblasts, as they transmit and amplify the traction forces generated by

The normalized solutions in time at **x** = **0**, furthest in the wound, are shown in Figure 6. This gives an indication of the development of the solutions in time. We see that the fibroblast concentration increases towards it equilibrium, whereas the myofibroblast concentration first rises and then falls again. This is due to the fact that, as the wound heals, the myofibroblasts either differentiate back to fibroblasts again or undergo apoptosis and hence eventually will disappear completely. Also the growth factor concentration eventually goes to zero as the wound is healed. The ECM density slowly rises in time and will eventually reach its equilibrium, allthough it can take some time before the ECM is completely restored. We finally

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12

x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

y

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

x

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

the fibroblasts, can clearly be seen.

Fig. 5. Displacement of the ECM two days after injury.

y

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

injury.

y

In this section, we present the model for angiogenesis. These two models each take a very different approach on how to model this process. Where the model due to Maggelakis (2003) focusses on the relation between a lack of oxygen and capillary growth, the model due to Gaffney (2002), attempts to model the migration of endothelial cells into the wound. Both models adress an important aspect of angiogenesis, but on the other hand both also miss an important aspect.

In this section we present a novel angiogenesis model based on the models of Maggelakis and Gaffney et al. In this model we attempt to unite the two approaches and thus create a model that is suitable for dealing with both aspects of angiogenesis. This gives a new model for angiogenesis. By combining these two aspects in one model we hope to create a more accurate model of the process of angiogenesis. As a basis we take the model of Gaffney et al. and extend it with a negative feedback mechanism for oxygen. This assures that both the endothelial cell migration and the oxygen shortage aspects are covered in this novel angiogenesis model.

The model due to Gaffney et al. is given by

$$\begin{aligned} \frac{\partial u\_{\text{tip}}}{\partial t} - \nabla \cdot (D\_1 \nabla u\_{\text{tip}} + D\_2 \nabla u\_{\text{end}}) &= f(u\_{\text{tip}}, u\_{\text{end}}), \\\\ \frac{\partial u\_{\text{end}}}{\partial t} - \nabla \cdot \lambda\_5 (D\_1 \nabla u\_{\text{tip}} + D\_2 \nabla u\_{\text{end}}) &= g(u\_{\text{tip}}, u\_{\text{end}}) \end{aligned} \tag{14}$$

Contraction and Angiogenesis 13

A Mathematical Model for Wound Contraction and Angiogenesis 501

λ2 λ6

<sup>2</sup> <sup>=</sup> 0.83 and *<sup>λ</sup>*<sup>0</sup>

end − *u*end)

**<sup>x</sup>**| |**x**| ≤ <sup>1</sup> 2 .

<sup>6</sup> = 1.

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

MDGF concentration



tip − *λ*4*u*tip*u*end,

end <sup>−</sup> *<sup>u</sup>*end) + *<sup>χ</sup>u*tip*u*end(*u*<sup>1</sup>

blood flow in the tips and capillaries (via diffusion through the capillary walls);

*<sup>u</sup>*tip <sup>−</sup> *<sup>λ</sup>*3*u*<sup>2</sup>

*au*end(*u*<sup>0</sup>

tip + *λ*4*u*tip*u*end).

To illustrate how this new model behaves we show some results in Figures 8 to 10. The input-data can be found in the appendix. As a computational domain we use the unit square,

In Figure 8 we see the oxygen and MDGF concentration seven days after injury. The oxygen tension decreases as one moves into the center of the wound. The relation between a lack of oxygen and the production of the macrophage derived growth factor can clearly be seen. In

i.e. Ω = {**x** = (*x*, *y*)| 0 ≤ *x*, *y* ≤ 1}, and the initial wound is given by Ω*<sup>w</sup>* =

Fig. 7. *λ*<sup>2</sup> and *λ*<sup>6</sup> as functions of *c*md. Here we used *τ*tip = *τ*end = 1, *λ*<sup>0</sup>

secretion by the macrophages as a result of a shortage of oxygen;

The last two equations were described earlier. Further, the functions *f* and *g* are given by

2

6

*c*md *c*md + *τ*tip

*c*md *c*md + *τ*end

+ *λ*5(*λ*3*u*<sup>2</sup>

*f*(*u*tip, *u*end) =*λ*<sup>0</sup>

*g*(*u*tip, *u*end) =*λ*<sup>0</sup>

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

where the left-hand sides account for accumulation and transport of tips and endothelial cells by means of a biased random walk process. Further, the functions *f* and *g* are given by

$$f(\boldsymbol{u}\_{\rm tip}, \boldsymbol{u}\_{\rm end}) = \lambda\_2 \boldsymbol{u}\_{\rm tip} - \lambda\_3 \boldsymbol{u}\_{\rm tip}^2 - \lambda\_4 \boldsymbol{u}\_{\rm tip} \boldsymbol{u}\_{\rm end},$$

$$g(\boldsymbol{u}\_{\rm tip}, \boldsymbol{u}\_{\rm end}) = \lambda\_6 a \boldsymbol{u}\_{\rm end} (\boldsymbol{u}\_{\rm end}^0 - \boldsymbol{u}\_{\rm end}) + \lambda\_6 \chi \boldsymbol{u}\_{\rm tip} \boldsymbol{u}\_{\rm end} (\boldsymbol{u}\_{\rm end}^1 - \boldsymbol{u}\_{\rm end}),$$

$$+ \lambda\_5 (\lambda\_3 \boldsymbol{u}\_{\rm tip}^2 + \lambda\_4 \boldsymbol{u}\_{\rm tip} \boldsymbol{u}\_{\rm end}).$$

These functions were taken from Gaffney (2002). The first term of the function *f* models tip-branching, as each tip has a likelihood to bifurcate. The second term of the function *f* accounts for tip-tip anastomosis, which is the process that two tips have a probability to join and hence form a loop, by which a tip is no longer a tip, and hence the number of tips decreases then. The third term of the function *f* stands for tip-sprout anastomosis, which is the process that a tip may converge into a branch, which closes the network as well and thereby also decreasing the number of tips. For the function *g*, we distinguish the following terms: The first one contains ordinary logistic proliferation under standard conditions, the second term accounts for an increased logistic growth due to the presence of tips. Finally, the third term accounts for the extent of anastomosis. To incorporate the effects of the macrophage derived growth factor (MDGF) concentration *c*md on the growth of capillary tips and the proliferation of endothelial cells we assume that *λ*<sup>2</sup> and *λ*<sup>6</sup> are functions of the macrophage derived growth factor *c*md, which is released as a result of a shortage of oxygen. These functions must be zero if there is no MDGF present and rise as the MDGF concentration rises. Furthermore the effects of additional MDGF must be lower if there is already a lot of MDGF present. Therefore we let

$$
\begin{aligned}
\lambda\_2(c\_{\rm md}) &= \lambda\_2^0 \frac{c\_{\rm md}}{c\_{\rm md} + \tau\_{\rm tip}}, \\
\lambda\_6(c\_{\rm md}) &= \lambda\_6^0 \frac{c\_{\rm md}}{c\_{\rm md} + \tau\_{\rm end}}.
\end{aligned}
$$

where *τ*tip and *τ*end qualify how respectively *λ*<sup>2</sup> and *λ*6, accounting for tip regeneration and logistic proliferation of endothelial cells respectively, depend on the MDGF concentration. The functions *λ*2(*c*md) and *λ*6(*c*md) are given in Figure 7 for 0 ≤ *c*md ≤ 1.

This settles the trigger mechanism that the MDGF concentration fullfills in the growth of new capillaries. To get a good overview of what the model looks like we give the PDEs that drive the model, i.e.

$$\frac{\partial \boldsymbol{u}\_{\rm{oxy}}}{\partial t} = D\_{\rm{oxy}} \Delta \boldsymbol{u}\_{\rm{oxy}} - \lambda\_{\rm{oxy}} \mu\_{\rm{oxy}} + \lambda\_{13} \left( \boldsymbol{u}\_{\rm{tip}} + \frac{\lambda\_{\rm{oxy}} \mu\_{\theta} \mu\_{\rm{end}}}{\lambda\_{13} \boldsymbol{u}\_{\rm{end}}^{0}} \right) \tag{15}$$

$$\frac{\partial c\_{\rm md}}{\partial t} = D\_{\rm md} \Delta c\_{\rm md} - \lambda\_{\rm md} c\_{\rm md} + \lambda\_{21} Q(u\_{\rm oxy}), \tag{16}$$

$$\frac{\partial \mu\_{\rm tip}}{\partial t} = \nabla \cdot \left\{ D\_1 \nabla \mu\_{\rm tip} + D\_2 \mu\_{\rm tip} \nabla \mu\_{\rm end} \right\} + f(\mu\_{\rm tip}, \mu\_{\rm end})\_\prime \tag{17}$$

$$\frac{\partial \mu\_{\rm end}}{\partial t} = \lambda\_1 \nabla \cdot \left\{ D\_1 \nabla u\_{\rm tip} + D\_2 u\_{\rm tip} \nabla u\_{\rm end} \right\} + g(u\_{\rm tip}, u\_{\rm end})\_\prime \tag{18}$$

where we explain the terms occuring in the above PDEs.

12 Will-be-set-by-IN-TECH

where the left-hand sides account for accumulation and transport of tips and endothelial cells by means of a biased random walk process. Further, the functions *f* and *g* are given by

tip − *λ*4*u*tip*u*end,

tip + *λ*4*u*tip*u*end).

These functions were taken from Gaffney (2002). The first term of the function *f* models tip-branching, as each tip has a likelihood to bifurcate. The second term of the function *f* accounts for tip-tip anastomosis, which is the process that two tips have a probability to join and hence form a loop, by which a tip is no longer a tip, and hence the number of tips decreases then. The third term of the function *f* stands for tip-sprout anastomosis, which is the process that a tip may converge into a branch, which closes the network as well and thereby also decreasing the number of tips. For the function *g*, we distinguish the following terms: The first one contains ordinary logistic proliferation under standard conditions, the second term accounts for an increased logistic growth due to the presence of tips. Finally, the third term accounts for the extent of anastomosis. To incorporate the effects of the macrophage derived growth factor (MDGF) concentration *c*md on the growth of capillary tips and the proliferation of endothelial cells we assume that *λ*<sup>2</sup> and *λ*<sup>6</sup> are functions of the macrophage derived growth factor *c*md, which is released as a result of a shortage of oxygen. These functions must be zero if there is no MDGF present and rise as the MDGF concentration rises. Furthermore the effects of additional MDGF must be lower if there is already a lot of MDGF present. Therefore we let

*<sup>λ</sup>*2(*c*md) = *<sup>λ</sup>*<sup>0</sup>

*λ*6(*c*md) = *λ*<sup>0</sup>

functions *λ*2(*c*md) and *λ*6(*c*md) are given in Figure 7 for 0 ≤ *c*md ≤ 1.

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>D</sup>*oxyΔ*u*oxy <sup>−</sup> *<sup>λ</sup>*oxy*u*oxy <sup>+</sup> *<sup>λ</sup>*<sup>13</sup>

*D*1∇*u*tip + *D*2*u*tip∇*u*end

*D*1∇*u*tip + *D*2*u*tip∇*u*end

the model, i.e.

*∂u*oxy

*∂c*md

*∂u*tip

*∂u*end

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> ∇ ·

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>λ</sup>*1∇ ·

where we explain the terms occuring in the above PDEs.

2

6

where *τ*tip and *τ*end qualify how respectively *λ*<sup>2</sup> and *λ*6, accounting for tip regeneration and logistic proliferation of endothelial cells respectively, depend on the MDGF concentration. The

This settles the trigger mechanism that the MDGF concentration fullfills in the growth of new capillaries. To get a good overview of what the model looks like we give the PDEs that drive

*c*md *c*md + *τ*tip

*c*md *c*md + *τ*end

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>D</sup>*mdΔ*c*md <sup>−</sup> *<sup>λ</sup>*md*c*md <sup>+</sup> *<sup>λ</sup>*21*Q*(*u*oxy), (16)

*u*tip +

*λ*oxy*uθu*end *λ*13*u*<sup>0</sup> end

+ *<sup>f</sup>*(*u*tip, *<sup>u</sup>*end), (17)

+ *g*(*u*tip, *u*end), (18)

, (15)

,

,

end <sup>−</sup> *<sup>u</sup>*end) + *<sup>λ</sup>*6*χu*tip*u*end(*u*<sup>1</sup>

end − *u*end)

*<sup>f</sup>*(*u*tip, *<sup>u</sup>*end) =*λ*2*u*tip <sup>−</sup> *<sup>λ</sup>*3*u*<sup>2</sup>

+ *λ*5(*λ*3*u*<sup>2</sup>

*g*(*u*tip, *u*end) =*λ*6*au*end(*u*<sup>0</sup>

Fig. 7. *λ*<sup>2</sup> and *λ*<sup>6</sup> as functions of *c*md. Here we used *τ*tip = *τ*end = 1, *λ*<sup>0</sup> <sup>2</sup> <sup>=</sup> 0.83 and *<sup>λ</sup>*<sup>0</sup> <sup>6</sup> = 1.


The last two equations were described earlier.

Further, the functions *f* and *g* are given by

$$\begin{split} f(\boldsymbol{u}\_{\text{tip}},\boldsymbol{u}\_{\text{end}}) &= \lambda\_{2}^{0}\frac{\boldsymbol{c}\_{\text{md}}}{\boldsymbol{c}\_{\text{md}}+\boldsymbol{\tau}\_{\text{tip}}}\boldsymbol{u}\_{\text{tip}} - \lambda\_{3}\boldsymbol{u}\_{\text{tip}}^{2} - \lambda\_{4}\boldsymbol{u}\_{\text{tip}}\boldsymbol{u}\_{\text{end}} \\ \boldsymbol{g}(\boldsymbol{u}\_{\text{tip}},\boldsymbol{u}\_{\text{end}}) &= \lambda\_{6}^{0}\frac{\boldsymbol{c}\_{\text{md}}}{\boldsymbol{c}\_{\text{md}}+\boldsymbol{\tau}\_{\text{end}}} \left(\boldsymbol{a}\boldsymbol{u}\_{\text{end}}(\boldsymbol{u}\_{\text{end}}^{0}-\boldsymbol{u}\_{\text{end}}) + \boldsymbol{\chi}\boldsymbol{u}\_{\text{tip}}\boldsymbol{u}\_{\text{end}}(\boldsymbol{u}\_{\text{end}}^{1}-\boldsymbol{u}\_{\text{end}})\right) \\ &+ \lambda\_{5}(\lambda\_{3}\boldsymbol{u}\_{\text{tip}}^{2}+\lambda\_{4}\boldsymbol{u}\_{\text{tip}}\boldsymbol{u}\_{\text{end}}). \end{split}$$

To illustrate how this new model behaves we show some results in Figures 8 to 10. The input-data can be found in the appendix. As a computational domain we use the unit square, i.e. Ω = {**x** = (*x*, *y*)| 0 ≤ *x*, *y* ≤ 1}, and the initial wound is given by Ω*<sup>w</sup>* = **<sup>x</sup>**| |**x**| ≤ <sup>1</sup> 2 .

In Figure 8 we see the oxygen and MDGF concentration seven days after injury. The oxygen tension decreases as one moves into the center of the wound. The relation between a lack of oxygen and the production of the macrophage derived growth factor can clearly be seen. In

x

Contraction and Angiogenesis 15

A Mathematical Model for Wound Contraction and Angiogenesis 503

Oxygen MDGF Capillary tips Endothelial cells

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> <sup>100</sup> <sup>0</sup>

Time in days

For all the stages during the proliferative phase of wound healing mathematical models have been presented that attempt to describe the processes involved in these stages. These models only focus on one wound healing stage particularly and do not take into account the interactions between these stages. One example of interaction is the need of oxygen for the growth of cells and thus the dependency of the profileration terms on the oxygen concentration. But there are several more interactions that we can think of, including some

In Vermolen (2010) and Vermolen (2011) an attempt is made to combine models of the three wound healing stages into one mathematical model. In this chapter we take the novel model on angiogenesis, presented the previous subsection, and make an attempt to couple it with the wound contraction model due to Javierre (2009), given in Section 2.1. We will investigate both the influence of the oxygen concentration on the wound contraction mechanism and the influence of the fibroblast concentration on capillary growth. Furthermore the displacement of the extracellular matrix (ECM) will also effect the angiogenesis process, this will also be

In the wound contraction stage several processes require energy to work. The growth of new tissue consumes energy as cells divide and form new cells, i.e. mitosis. Also the active movement of the fibroblasts requires energy since the cells crawl over each other. This energy needed for wound contraction comes primarly from the food that an individual consumes, however, oxygen is needed to put the energy to use, i.e. the cells require oxygen to get the most out of the energy. This means that the oxygen concentration is good indicator for the

Fig. 10. Normalized solutions furthest in the wound, at **x** = **0**, in time.

taken into account. For now we will not consider wound closure.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

that will not be taken into account here either.

amount of energy that is available.

**3. The coupled model**

Fig. 8. Normalized oxygen (left) and macrophage derived growth factor (right) concentration seven days after injury.

x

Fig. 9. Normalized capillary tip concentration (left) and endothelial cell density (right) seven days after injury.

areas where the oxygen concentration is low, i.e. inside the wound, the MDGF concentration is at its peak and vice versa. Hence, the MDGF concentration is maximal at the center of the wound. The time for ingress of macrophages has been assumed to be negligible. Also in Figure 10 we see that as the oxygen concentration rises the MDGF concentration drops at approximately the same rate.

Furthermore in Figure 9 we see a front of capillary tips, which slowly moves towards the center of the wound. Also the endothelial cell density is slighty elevated at the wound edge, since this is where the new capillaries are formed. This can also be seen in Figure 10, where the endothelial cell density first peaks and then drops again towards an equilibrium. The capillary tip concentration also peaks, but then drops to zero again. This is to be expected since eventually all capillary tips join together in the newly formed capillary network. The local maximum of the oxygen concentration at *t* ≈ 18 days, follows from the peak in capillary tips. As the shortage of oxygen decreases, the concentration of macrophage derived growth factors decreases down to zero.

Fig. 10. Normalized solutions furthest in the wound, at **x** = **0**, in time.

#### **3. The coupled model**

14 Will-be-set-by-IN-TECH

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

0.05 0.1 0.15 0.2 0.25 0.3

Fig. 9. Normalized capillary tip concentration (left) and endothelial cell density (right) seven

areas where the oxygen concentration is low, i.e. inside the wound, the MDGF concentration is at its peak and vice versa. Hence, the MDGF concentration is maximal at the center of the wound. The time for ingress of macrophages has been assumed to be negligible. Also in Figure 10 we see that as the oxygen concentration rises the MDGF concentration drops at

Furthermore in Figure 9 we see a front of capillary tips, which slowly moves towards the center of the wound. Also the endothelial cell density is slighty elevated at the wound edge, since this is where the new capillaries are formed. This can also be seen in Figure 10, where the endothelial cell density first peaks and then drops again towards an equilibrium. The capillary tip concentration also peaks, but then drops to zero again. This is to be expected since eventually all capillary tips join together in the newly formed capillary network. The local maximum of the oxygen concentration at *t* ≈ 18 days, follows from the peak in capillary tips. As the shortage of oxygen decreases, the concentration of macrophage derived growth

Fig. 8. Normalized oxygen (left) and macrophage derived growth factor (right) concentration

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.1 0.2 0.3 0.4 0.5 0.6 0.7

x

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

y

y

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

x

<sup>0</sup> 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 <sup>1</sup> <sup>0</sup>

x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

seven days after injury.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

days after injury.

approximately the same rate.

factors decreases down to zero.

y

y

For all the stages during the proliferative phase of wound healing mathematical models have been presented that attempt to describe the processes involved in these stages. These models only focus on one wound healing stage particularly and do not take into account the interactions between these stages. One example of interaction is the need of oxygen for the growth of cells and thus the dependency of the profileration terms on the oxygen concentration. But there are several more interactions that we can think of, including some that will not be taken into account here either.

In Vermolen (2010) and Vermolen (2011) an attempt is made to combine models of the three wound healing stages into one mathematical model. In this chapter we take the novel model on angiogenesis, presented the previous subsection, and make an attempt to couple it with the wound contraction model due to Javierre (2009), given in Section 2.1. We will investigate both the influence of the oxygen concentration on the wound contraction mechanism and the influence of the fibroblast concentration on capillary growth. Furthermore the displacement of the extracellular matrix (ECM) will also effect the angiogenesis process, this will also be taken into account. For now we will not consider wound closure.

In the wound contraction stage several processes require energy to work. The growth of new tissue consumes energy as cells divide and form new cells, i.e. mitosis. Also the active movement of the fibroblasts requires energy since the cells crawl over each other. This energy needed for wound contraction comes primarly from the food that an individual consumes, however, oxygen is needed to put the energy to use, i.e. the cells require oxygen to get the most out of the energy. This means that the oxygen concentration is good indicator for the amount of energy that is available.

Contraction and Angiogenesis 17

A Mathematical Model for Wound Contraction and Angiogenesis 505

respectively. Similarly to how the oxygen concentration is coupled with the fibroblast PDE we couple the fibroblasts to the capillary tip and endothelial cell density PDEs. Therefor we

> *u*fib/*u*<sup>0</sup> fib *τ*fib + *u*fib/*u*<sup>0</sup>

in front of both production terms. This results in no capillary growth if there are no fibroblasts present, i.e. there is no appropriate dermal tissue for them to grow in. As the fibroblast concentration rises the rate of capillary growth also rises towards its normal rate, since the

This concludes the coupling between the wound contraction model due to Javierre and the novel angiogenesis model covered in this thesis. Of course there are several other interaction between them. One can for instance think about the differentation from fibroblasts to myofibroblasts and back, which surely also consumes energy. Alternatively, we could think of also linking the chemotaxis term to the oxygen content. This could be a subject for further

Next we present some results of the computation done on the coupled model. The values of the parameters used can be found in appendix. We vary the diffusion speeds of the fibroblasts and the oxygen, since we consider these to be of great importance (the problem seems to be diffusion dominated) and we would like to study the effect of these parameters on the

In Figures 11 to 13 the results of the simulations can be found. The figures show the normalized solutions in time furthest in the wound, i.e. at **x** = **0**. The difference in diffusion

We see that with slow fibroblast diffusion the growth of capillaries start off at a later point in time. This is to be expected since the fibroblasts provide the tissue in which the capillaries can grow. The growth of the capillaries does follow a similar course. First the capillary tips find their way to the center of the wound and after that the capillary network is restored, which decreases the number of tips. As the number of tips decreases down to zero, the oxygen content decreases for a short while due to the lack of this tip-source. Later, the oxygen tension

Further, having slower oxygen diffusion, the oxygen concentration depends more on the growth of new capillaries. In Figure 13 we see that the oxygen concentration rises far slower than 11 and 12. The oxygen concentration grows due to capillaries transporting oxygen to the

We see that the varying the diffusion speeds has a major effect on the solutions. This confirms our idea that the problem is diffusion dominated. Although the solutions follow similar

Furthermore the coupling between the two models can also be seen. Especially with low diffusion speeds we see that the fibroblast concentration clearly depends on the oxygen concentration. Also the growth of new capillaries starts off later than in the uncoupled model, see Section 2.2. This can also be explained by the coupling, since the growth of capillaries now

end <sup>−</sup> *<sup>u</sup>*end) + *<sup>χ</sup>u*tip*u*end(*u*<sup>1</sup>

fib

end − *u*end)

and

study.

solution.

speeds can clearly be seen in the graphs.

increases to its undamaged equilibrium.

wound side and not primarly due to diffusion.

patterns in all cases, the differences can also be seen clearly.

*λ*0 6

introduce the factor

term approaches one.

*c*md *c*md + *τ*end *au*end(*u*<sup>0</sup>

If we look at the partial differential equation (PDE) for the fibroblast concentration in the Javierre model,

$$\frac{\partial u\_{\rm fib}}{\partial t} + \nabla \cdot \left( \frac{\partial \mathbf{u}}{\partial t} u\_{\rm fib} - D\_{\rm fib} \nabla u\_{\rm fib} + \frac{a\_{\rm fib}}{(b\_{\rm fib} + c\_{\rm ferm})^2} u\_{\rm fib} \nabla c\_{\rm ferm} \right) =$$

$$\left( \lambda\_{\rm fib} + \frac{\lambda\_{\rm fib}^0 c\_{\rm ferm}}{\mathbb{C}\_{1/2} + c\_{\rm ferm}} \right) u\_{\rm fib} \left( 1 - \frac{u\_{\rm fib}}{K} \right) - \frac{p\_{\rm cell}(\theta)}{\tau\_d + p\_{\rm cell}(\theta)} \frac{k\_1 c\_{\rm ferm}}{\mathbb{C}\_k + c\_{\rm ferm}} u\_{\rm fib} + k\_2 u\_{\rm mov} - d\_{\rm fib} u\_{\rm fib} \tag{19}$$

the third term on the left hand side denotes the active mobility due to random walk and the first term on the right hand side represents fibroblast proliferation. To incorporate the effects of the oxygen level on these two processes we replace these terms with the following,

$$\frac{\mu\_{\rm{oxy}}}{\pi\_{\rm{oxy}} + \mu\_{\rm{oxy}}} (-D\_{\rm{fib}} \nabla \mu\_{\rm{fib}} + \frac{a\_{\rm{fib}}}{(b\_{\rm{fib}} + c\_{\rm{ecm}})^2} \mu\_{\rm{fib}} \nabla c\_{\rm{ecm}})^2$$

for the mobility and

$$\frac{\mu\_{\rm{oxy}}}{\pi\_{\rm{oxy}} + \mu\_{\rm{oxy}}} \left(\lambda\_{\rm{fib}} + \frac{\lambda\_{\rm{fib}}^{0} c\_{\rm{ecm}}}{C\_{1/2} + c\_{\rm{ecm}}}\right) \mu\_{\rm{fib}} \left(1 - \frac{\mu\_{\rm{fib}}}{K}\right).$$

for the proliferation of fibroblasts. The extra term in front, *<sup>u</sup>*oxy *<sup>τ</sup>*oxy+*u*oxy , assures that if there is no oxygen present, i.e. no energy can be consumed, the two processes are stopped. When the tissue is saturated with oxygen the processes continue at their normal rate (the term then approaches one). Note that due to the PDE for the oxygen concentration the oxygen level will never rise to dangerous levels (the concentration moves towards an equilibrium), i.e. oxygen poisoning will never take place.

The proliferation of myofibroblasts also consumes energy and thus the same dependency on the oxygen concentration can be found in the PDE for the myofibroblasts. This concludes the influence of the oxygen level on wound contraction covered in this model.

The other way around wound contraction also effects the growth of new capillaries. First the displacement of the ECM causes passive convection of the variables of the angiogenesis model. This means that in each of the four equations an extra term is found that describes this passive convection, i.e.

$$\nabla \cdot \left(\frac{\partial \mathbf{u}}{\partial t} u\_i\right)\_{\prime \prime}$$

where *ui* denote the four variables of the novel angiogenesis model.

Furthermore the new capillaries need tissue to grow in. At first all tissue in the wound has been destroyed by the injury and so the capillaries can not grow. Gradually fibroblasts invade the wound and new tissue is formed. Only then can the recovery of the capillary network start. This is why it is reasonable to let the growth of capillaries depend on the fibroblast concentration.

The growth terms in (17) and (18) are given by

$$\lambda\_2^0 \frac{c\_{\rm md}}{c\_{\rm md} + \tau\_{\rm tip}} u\_{\rm tip}$$

and

16 Will-be-set-by-IN-TECH

If we look at the partial differential equation (PDE) for the fibroblast concentration in the

<sup>−</sup> *<sup>p</sup>*cell(*θ*) *τ<sup>d</sup>* + *p*cell(*θ*)

the third term on the left hand side denotes the active mobility due to random walk and the first term on the right hand side represents fibroblast proliferation. To incorporate the effects

> *λ*0 fib*c*ecm *C*1/2 + *c*ecm

no oxygen present, i.e. no energy can be consumed, the two processes are stopped. When the tissue is saturated with oxygen the processes continue at their normal rate (the term then approaches one). Note that due to the PDE for the oxygen concentration the oxygen level will never rise to dangerous levels (the concentration moves towards an equilibrium), i.e. oxygen

The proliferation of myofibroblasts also consumes energy and thus the same dependency on the oxygen concentration can be found in the PDE for the myofibroblasts. This concludes the

The other way around wound contraction also effects the growth of new capillaries. First the displacement of the ECM causes passive convection of the variables of the angiogenesis model. This means that in each of the four equations an extra term is found that describes this

> *∂***u** *∂t ui* ,

Furthermore the new capillaries need tissue to grow in. At first all tissue in the wound has been destroyed by the injury and so the capillaries can not grow. Gradually fibroblasts invade the wound and new tissue is formed. Only then can the recovery of the capillary network start. This is why it is reasonable to let the growth of capillaries depend on the fibroblast

> *c*md *c*md + *τ*tip

*u*tip

∇ ·

*λ*0 2

of the oxygen level on these two processes we replace these terms with the following,

(−*D*fib∇*u*fib +

*λ*fib +

influence of the oxygen level on wound contraction covered in this model.

where *ui* denote the four variables of the novel angiogenesis model.

The growth terms in (17) and (18) are given by

for the proliferation of fibroblasts. The extra term in front, *<sup>u</sup>*oxy

*u*fib − *D*fib∇*u*fib +

*a*fib

*k*1*c*ecm *Ck* + *c*ecm

(*b*fib <sup>+</sup> *<sup>c</sup>*ecm)<sup>2</sup> *<sup>u</sup>*fib∇*c*ecm)

*a*fib

 *u*fib <sup>1</sup> <sup>−</sup> *<sup>u</sup>*fib *K* 

(*b*fib <sup>+</sup> *<sup>c</sup>*ecm)<sup>2</sup> *<sup>u</sup>*fib∇*c*ecm

 =

*u*fib + *k*2*u*myo − *d*fib*u*fib, (19)

*<sup>τ</sup>*oxy+*u*oxy , assures that if there is

 *∂***u** *∂t*

Javierre model,

*λ*fib +

*λ*0 fib*c*ecm *C*1/2 + *c*ecm

for the mobility and

poisoning will never take place.

passive convection, i.e.

concentration.

*∂u*fib *<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ ·

*u*oxy *τ*oxy + *u*oxy

*u*oxy *τ*oxy + *u*oxy

 *u*fib <sup>1</sup> <sup>−</sup> *<sup>u</sup>*fib *K* 

$$
\lambda\_6^0 \frac{c\_{\rm md}}{c\_{\rm md} + \tau\_{\rm end}} \left( a u\_{\rm end} (\mu\_{\rm end}^0 - u\_{\rm end}) + \chi u\_{\rm tip} \mu\_{\rm end} (\mu\_{\rm end}^1 - u\_{\rm end}) \right),
$$

respectively. Similarly to how the oxygen concentration is coupled with the fibroblast PDE we couple the fibroblasts to the capillary tip and endothelial cell density PDEs. Therefor we introduce the factor

$$\frac{\mu\_{\rm fib}/\mu\_{\rm fib}^0}{\tau\_{\rm fib} + \mu\_{\rm fib}/\mu\_{\rm fib}^0}$$

in front of both production terms. This results in no capillary growth if there are no fibroblasts present, i.e. there is no appropriate dermal tissue for them to grow in. As the fibroblast concentration rises the rate of capillary growth also rises towards its normal rate, since the term approaches one.

This concludes the coupling between the wound contraction model due to Javierre and the novel angiogenesis model covered in this thesis. Of course there are several other interaction between them. One can for instance think about the differentation from fibroblasts to myofibroblasts and back, which surely also consumes energy. Alternatively, we could think of also linking the chemotaxis term to the oxygen content. This could be a subject for further study.

Next we present some results of the computation done on the coupled model. The values of the parameters used can be found in appendix. We vary the diffusion speeds of the fibroblasts and the oxygen, since we consider these to be of great importance (the problem seems to be diffusion dominated) and we would like to study the effect of these parameters on the solution.

In Figures 11 to 13 the results of the simulations can be found. The figures show the normalized solutions in time furthest in the wound, i.e. at **x** = **0**. The difference in diffusion speeds can clearly be seen in the graphs.

We see that with slow fibroblast diffusion the growth of capillaries start off at a later point in time. This is to be expected since the fibroblasts provide the tissue in which the capillaries can grow. The growth of the capillaries does follow a similar course. First the capillary tips find their way to the center of the wound and after that the capillary network is restored, which decreases the number of tips. As the number of tips decreases down to zero, the oxygen content decreases for a short while due to the lack of this tip-source. Later, the oxygen tension increases to its undamaged equilibrium.

Further, having slower oxygen diffusion, the oxygen concentration depends more on the growth of new capillaries. In Figure 13 we see that the oxygen concentration rises far slower than 11 and 12. The oxygen concentration grows due to capillaries transporting oxygen to the wound side and not primarly due to diffusion.

We see that the varying the diffusion speeds has a major effect on the solutions. This confirms our idea that the problem is diffusion dominated. Although the solutions follow similar patterns in all cases, the differences can also be seen clearly.

Furthermore the coupling between the two models can also be seen. Especially with low diffusion speeds we see that the fibroblast concentration clearly depends on the oxygen concentration. Also the growth of new capillaries starts off later than in the uncoupled model, see Section 2.2. This can also be explained by the coupling, since the growth of capillaries now

Contraction and Angiogenesis 19

A Mathematical Model for Wound Contraction and Angiogenesis 507

Oxygen Fibroblasts Capillary tips Endothelial cells

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> <sup>100</sup> −0.1

depends on the fibroblast concentration. So the growth can not begin untill there is enough

We use a Galerkin finite-element method with linear triangular elements to solve the resulting system of PDEs. The reaction-transport equations are solved using a second-order Trapezoidal Time Numerical Integration Method. The nonlinear problem at each time step is solved using a Newton's method, which gives quadratic convergence. For the numerical integration of the spatial integrals in the weak formulation, we use the second order accurate Newton-Cotes integration rule. The numerical solution will thereby have a quadratic accuracy in both element size and time step. For an efficient implementation, the reactive-transport equations are solved as a fully coupled problem in which the degrees of freedom are numbered in a nodewise fashion. The mechanical parameters were taken at the current time step, since at each time-step the mechanical problem is solved before the biological reactive-transport

The visco-elastic equations are solved using a Galerkin finite-element method with linear triangular elements. The system is solved as a fully coupled problem with a nodewise arrangement of the degrees of freedom, being the vertical and horizontal displacement. The concentrations and cellular densities from the previous time step are used in the mechanical balance equations from visco-elasticity. For the time-integration, a Trapezoidal Numerical Integration Method is used as well to have second order accuracy in time as well. The linear

Fig. 13. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.002

Time in days

0

cm2/days and *D*oxy = 0.001 cm2/days.

**4. Numerical methods**

equations are solved numerically.

tissue available, i.e. the fibroblast concentration is non-zero.

elements give a second-order accuracy for the element size.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Fig. 11. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.02 cm2/days and *D*oxy = 0.01 cm2/days.

Fig. 12. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.002 cm2/days and *D*oxy = 0.01 cm2/days.

18 Will-be-set-by-IN-TECH

Oxygen Fibroblasts Capillary tips Endothelial cells

Oxygen Fibroblasts Capillary tips Endothelial cells

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> <sup>100</sup> <sup>0</sup>

Time in days

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> <sup>100</sup> <sup>0</sup>

Time in days

Fig. 12. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.002

Fig. 11. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.02

0.1

cm2/days and *D*oxy = 0.01 cm2/days.

0.1

cm2/days and *D*oxy = 0.01 cm2/days.

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 13. Normalized solutions furthest in the wound, at **x** = **0**, in time for *D*fib = 0.002 cm2/days and *D*oxy = 0.001 cm2/days.

depends on the fibroblast concentration. So the growth can not begin untill there is enough tissue available, i.e. the fibroblast concentration is non-zero.

#### **4. Numerical methods**

We use a Galerkin finite-element method with linear triangular elements to solve the resulting system of PDEs. The reaction-transport equations are solved using a second-order Trapezoidal Time Numerical Integration Method. The nonlinear problem at each time step is solved using a Newton's method, which gives quadratic convergence. For the numerical integration of the spatial integrals in the weak formulation, we use the second order accurate Newton-Cotes integration rule. The numerical solution will thereby have a quadratic accuracy in both element size and time step. For an efficient implementation, the reactive-transport equations are solved as a fully coupled problem in which the degrees of freedom are numbered in a nodewise fashion. The mechanical parameters were taken at the current time step, since at each time-step the mechanical problem is solved before the biological reactive-transport equations are solved numerically.

The visco-elastic equations are solved using a Galerkin finite-element method with linear triangular elements. The system is solved as a fully coupled problem with a nodewise arrangement of the degrees of freedom, being the vertical and horizontal displacement. The concentrations and cellular densities from the previous time step are used in the mechanical balance equations from visco-elasticity. For the time-integration, a Trapezoidal Numerical Integration Method is used as well to have second order accuracy in time as well. The linear elements give a second-order accuracy for the element size.

Contraction and Angiogenesis 21

A Mathematical Model for Wound Contraction and Angiogenesis 509

The first important innovations in the present study are the formulation of a new angiogenesis model, which consists of a driving force due to lack of oxygen (Maggelakis) and a balance of both the capillary tips and capillaries by tracking these quantities as individual parameters. The second important innovation concerns the combination of the wound contraction model due to Javierre (2009) with the newly developed angiogenesis model. As an important factor of communication between both processes, we use the oxygen tension (on fibroblast mobility and proliferation), the presence of macrophage derived growth factor influencing the logistic proliferation rate of endothelial cells and the fibroblast density on the production of endothelial cells. Of course, these coupling terms have been assumed from educated guesses. An experimental validation could be decisive on the validity of the assumptions used in the present model. The calculations in the present manuscript are preliminary and a parameter

We also note that there is a large uncertainty for the values of all the parameters used. We intend to use stochastic methods to better deal with the occurrance of patient dependencies and other uncertainties in the data. This will be tackled in future studies. The present paper

In this Appendix, we show the data that we used for the calculations. The first two tables show the parameter values for the skin regeneration models ((myo)fibroblasts and collagen). The third table gives the data of the visco-elastic model for wound contraction. The fourth table gives the input for the angiogenesis model. Further, the values of the coupling parameters are

Parameter Description Value Dimension

fib Fibroblast density in healthy tissue 104 cells/cm<sup>3</sup>

fib Fibroblast diffusion rate <sup>2</sup> <sup>×</sup> <sup>10</sup>−2/2 <sup>×</sup> <sup>10</sup>−<sup>3</sup> cm2/day *<sup>a</sup>*fib Determines, with *<sup>b</sup>*fib, the maximal chemotaxis rate per unit GF conc. 4 <sup>×</sup> <sup>10</sup>−<sup>10</sup> g/cm day *<sup>b</sup>*fib GF conc. that produces 25% of the maximal chemotatic response <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>9</sup> g/cm<sup>3</sup> *λ*fib Fibroblast proliferation rate 0.832 day−<sup>1</sup>

fib Maximal GF induced proliferation rate 0.3 day−<sup>1</sup> *<sup>K</sup>* Determines the fibroblast equilibrium density <sup>1</sup> <sup>×</sup> 107 cells/cm<sup>3</sup> *<sup>C</sup>*1/2 Half-maximal GF enhancement of fibroblasts proliferation <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>8</sup> g/cm<sup>3</sup> *k*<sup>1</sup> Maximal fibroblast to myofibroblast differentation rate 0.8 day−<sup>1</sup> *k*<sup>2</sup> Myofibroblast to fibroblast differentation rate 0.693 day−<sup>1</sup> *Ck* Half-maximal GF enhancement of fibroblast to myofibroblast differentation 10−<sup>8</sup> g/cm<sup>3</sup> *d*fib Fibroblast death rate 0.831 day−<sup>1</sup> *ε*myo Myofibroblast to fibroblast logistic growth rate proportionality factor 0.5 *<sup>d</sup>*myo Myofibroblasts death rate 2.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> day−<sup>1</sup>

*<sup>τ</sup>*oxy Determines oxygen enhancement of (myo)fibroblasts proliferation 2 <sup>×</sup> <sup>10</sup>−<sup>4</sup> Table 1. Parameters related to the fibroblast and myofibroblast equations.

sensitivity analysis is to be carried out.

**7. Appendix: Parameter values**

scattered throughout the tables as well.

*u*0

*D*<sup>0</sup>

*λ*0

was more descriptive, rather than mathematical.

#### **5. Discussion**

In this paper, we described and coupled some mathematical models for wound contraction and angiogenesis as subprocesses for dermal regeneration post-wounding. It can be seen in Figure 10 that in the course of healing, the endothelial cell density, which is the most important measure for the vascular network density, first increases up to value above the equilibrium value, which corresponds to the fully healed state. This increase coincides with the temporary high capillary tip density at these times. After this intermediate stage, the capillary tip density decreases down to zero, and the endothelial cell density decreases to the equilibrium as time proceeds. This qualitatively agrees with the rabbit ear chamber images from the experiments that were carried out by Komori (2005), among others.

Furthermore, the numerical solutions converge to the undamaged state at the longest scales. This is consistent with clinical situations if the initial wound is very small. However, scar formation due to an excessive regeneration of collagen is not taken into account in the present model. Hence, the present model is incapable of predicting the occurrance of ulcers or hypertrophic scars that can result for large wounds or burns. This issue will be dealt with in future, as well as the extension of the present model by the incorporation of the epidermal layer which starts healing as a result of triggering of keratinocytes by the signaling agents that are secreted by the fibroblasts in the dermis. We are also interested in the development of the basal membrane, which separates the dermis from epidermis, in the course of time.

We also note that the model is entirely based on the continuum hypothesis, and that all densities and concentrations should be considered as averaged quantities over a *representative volume element*. At the smaller scale, one has to use a more cellular approach, such as the cellular automata (cellular Potts) models or (semi-) stochastic approaches. A translation between these classes of models is highly desirable and hardly found in literature.

Finally, we note that modeling of wound healing is very complicated, due to the enormous complexity of the healing process. The present model can already help in obtaining insight into the relation between several parameters, such as the influence of various parameters on the speed and quality of wound healing. In future studies, we intend to include the effects of several treatment of wounds, such as pressure or drug treatments, on the kinetics of wound healing. The models can then be used as a tool for medical doctors to evaluate the influence of certain treatments on wound healing, as well as to determine the circumstances for optimal healing with respect to minimization of hypertrophic scars in the case of burns. Here, we also remark that all data are patient-dependent and hence all parameters in the model are exposed to a stochastic nature. In order to deal with the issues of uncertainties, we intend to employ stochastic finite element methods in near future.

#### **6. Conclusions**

We presented an advanced model for the coupling of angiogenesis and wound contraction. The model gives a more complete picture than most of the presently existing formalisms in literature. However, still some room for improvement remains. The current model could also be extended with the healing of the epidermis located on top of the dermis. This will be done in a future study.

20 Will-be-set-by-IN-TECH

In this paper, we described and coupled some mathematical models for wound contraction and angiogenesis as subprocesses for dermal regeneration post-wounding. It can be seen in Figure 10 that in the course of healing, the endothelial cell density, which is the most important measure for the vascular network density, first increases up to value above the equilibrium value, which corresponds to the fully healed state. This increase coincides with the temporary high capillary tip density at these times. After this intermediate stage, the capillary tip density decreases down to zero, and the endothelial cell density decreases to the equilibrium as time proceeds. This qualitatively agrees with the rabbit ear chamber images from the experiments

Furthermore, the numerical solutions converge to the undamaged state at the longest scales. This is consistent with clinical situations if the initial wound is very small. However, scar formation due to an excessive regeneration of collagen is not taken into account in the present model. Hence, the present model is incapable of predicting the occurrance of ulcers or hypertrophic scars that can result for large wounds or burns. This issue will be dealt with in future, as well as the extension of the present model by the incorporation of the epidermal layer which starts healing as a result of triggering of keratinocytes by the signaling agents that are secreted by the fibroblasts in the dermis. We are also interested in the development of the

We also note that the model is entirely based on the continuum hypothesis, and that all densities and concentrations should be considered as averaged quantities over a *representative volume element*. At the smaller scale, one has to use a more cellular approach, such as the cellular automata (cellular Potts) models or (semi-) stochastic approaches. A translation

Finally, we note that modeling of wound healing is very complicated, due to the enormous complexity of the healing process. The present model can already help in obtaining insight into the relation between several parameters, such as the influence of various parameters on the speed and quality of wound healing. In future studies, we intend to include the effects of several treatment of wounds, such as pressure or drug treatments, on the kinetics of wound healing. The models can then be used as a tool for medical doctors to evaluate the influence of certain treatments on wound healing, as well as to determine the circumstances for optimal healing with respect to minimization of hypertrophic scars in the case of burns. Here, we also remark that all data are patient-dependent and hence all parameters in the model are exposed to a stochastic nature. In order to deal with the issues of uncertainties, we intend to employ

We presented an advanced model for the coupling of angiogenesis and wound contraction. The model gives a more complete picture than most of the presently existing formalisms in literature. However, still some room for improvement remains. The current model could also be extended with the healing of the epidermis located on top of the dermis. This will be done

basal membrane, which separates the dermis from epidermis, in the course of time.

between these classes of models is highly desirable and hardly found in literature.

that were carried out by Komori (2005), among others.

stochastic finite element methods in near future.

**6. Conclusions**

in a future study.

**5. Discussion**

The first important innovations in the present study are the formulation of a new angiogenesis model, which consists of a driving force due to lack of oxygen (Maggelakis) and a balance of both the capillary tips and capillaries by tracking these quantities as individual parameters. The second important innovation concerns the combination of the wound contraction model due to Javierre (2009) with the newly developed angiogenesis model. As an important factor of communication between both processes, we use the oxygen tension (on fibroblast mobility and proliferation), the presence of macrophage derived growth factor influencing the logistic proliferation rate of endothelial cells and the fibroblast density on the production of endothelial cells. Of course, these coupling terms have been assumed from educated guesses. An experimental validation could be decisive on the validity of the assumptions used in the present model. The calculations in the present manuscript are preliminary and a parameter sensitivity analysis is to be carried out.

We also note that there is a large uncertainty for the values of all the parameters used. We intend to use stochastic methods to better deal with the occurrance of patient dependencies and other uncertainties in the data. This will be tackled in future studies. The present paper was more descriptive, rather than mathematical.

#### **7. Appendix: Parameter values**

In this Appendix, we show the data that we used for the calculations. The first two tables show the parameter values for the skin regeneration models ((myo)fibroblasts and collagen). The third table gives the data of the visco-elastic model for wound contraction. The fourth table gives the input for the angiogenesis model. Further, the values of the coupling parameters are scattered throughout the tables as well.


Table 1. Parameters related to the fibroblast and myofibroblast equations.

Contraction and Angiogenesis 23

A Mathematical Model for Wound Contraction and Angiogenesis 511

Parameter Description Value Dimension

oxy Oxygen concentration in healthy tissue 5 mg/cm<sup>3</sup> *D*oxy Oxygen diffusion rate 10−2/10−<sup>3</sup> day−<sup>1</sup> *<sup>λ</sup>*oxy Oxygen decay rate <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> day−<sup>1</sup> *λ*<sup>13</sup> Oxygen transport rate 1 day−<sup>1</sup> *u<sup>θ</sup>* Threshold value for macrophage derived GF 5 mg/cm<sup>3</sup> *D*md Macrophage derived GF diffusion rate 0.1 day−<sup>1</sup> *λ*md Macrophage derived GF decay rate 1 day−<sup>1</sup> *λ*<sup>21</sup> Macrophage derived GF production rate 10 day−<sup>1</sup>

tip Normalized initial capillary tip concentration in the small strip facing the wound 1 - *<sup>D</sup>*<sup>1</sup> Capillary tip diffusion rate 3.5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> cm2/day

end Normalized endothelial cell density in healthy tissue 1 - *<sup>D</sup>*<sup>2</sup> Endothelial cell diffusion rate 3.5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> cm2/day

<sup>2</sup> Capillary tip growth rate 0.83 day−<sup>1</sup> *τ*tip Half-maximal macropahge derived GF enhancement of capillary tip growth 1 *λ*<sup>3</sup> Rate at which two capillary tips meet 0.83 day−<sup>1</sup> *λ*<sup>4</sup> Rate at which a capillary tip meets another capillary 0.85 day−<sup>1</sup>

<sup>6</sup> Endothelial cell proliferation rate 1 day−<sup>1</sup> *τ*end Half-maximal macropahge derived GF enhancement of endothelial cell proliferation 1 -

Bloemen, M.C.T., van der Veer, W.M. Ulrich, M.M.W., van Zuijlen, F.B., Niessen, F.B.,

van der Veer, W.M., Bloemen, M.C.T., Ulrich, M.M.W., Molema, G., van Zuijlen, P.P.M.,

Dallon, J.C., Ehrlich, H.P.: *A review of fibroblast populated collagen lattices* Wound repair and

Dallon, J.C.: *Multiscale modeling of cellular systems in biology*, Current opinion in colloid and

urray, J.D.: *Mathematical Biology II: Spatial Models and Biomedical Applications*, New York:

Olsen, L., Sherratt, J.A. and Maini, P.K.: *A mechanochemical model for adult dermal wound closure*

Maggelakis, A.: *A mathematical model for tissue replacement during epidermal wound healing*,

Gaffney, E.A., Pugh, K. and Maini, P.K.: *Investigating a simple model for cutaneous wound healing*

Sherratt, J.A. and Murray, J.D.: *Mathematical analysis of a basic model for epidermal wound healing*,

Adam, J.A.: *A simplified model of wound healing (with particular reference to the critical size defect)*,

Applied Mathematical Modelling 2003; 27:189-196.

Journal of Theoretical Biology 1991; 29:389-404.

*angiogenesis*, Journal of Theoretical Biology 2002; 45:337-374.

Mathematical and Computer Modelling 1999; 30:23-32.

*and the permanence of contracted tissue displacement*, Journal of Theoretical Biology 1995;

Middelkoop, E.: *Prevention and curative management of hypertrophic scar formation*, to

Middelkoop, E., Niessen, F.B.: *Potential cellular and molecular causes of hypertrophic*

end Determines, with *a* and *χ*, the equilibrium endothelial cell density 10 *λ*<sup>5</sup> Capillary tip to capillary proportionality factor 0.25 *τ*fib Half-maximal fibroblast enhancement of capillary growth and endothelial cell proliferation 1 -

end, the equilibrium endothelial cell density 0.25 -

end, the equilibrium endothelial cell density 0.3 -

*u*0

*u*0

*u*0

*λ*0

*λ*0

*u*1

**8. References**

*a* Determines, with *χ* and *u*<sup>1</sup>

*χ* Determines, with *a* and *u*<sup>1</sup>

Table 4. Parameters related to angiogenesis.

appear in Burns, 2011.

Springer-Verlag; 2004.

177:113-128.

*scar formation*, Burns, 35, 15–29, 2009.

regeneration, 16, 472–479, 2008.

interface science, 15, 24–31, 2010.






Table 4. Parameters related to angiogenesis.

#### **8. References**

22 Will-be-set-by-IN-TECH

Parameter Description Value Dimension *ρ*<sup>0</sup> Collagen concentration in healthy tissue 0.1 g/cm<sup>3</sup> *ρ*ini Initial collagen concentration in the wound 10−<sup>3</sup> g/cm<sup>3</sup> *λρ* Collagen production rate 7.59 <sup>×</sup> <sup>10</sup>−<sup>10</sup> <sup>g</sup>3/cm6 cell day

*<sup>ρ</sup>* Maximal rate of GF induced collagen production 7.59 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>g</sup>3/cm6 cell day *C<sup>ρ</sup>* Half-maximal GF enhancement of collagen production 10−<sup>8</sup> g/cm<sup>3</sup> *R<sup>ρ</sup>* Half-maximal collagen enhancement of ECM deposition 0.3 g/cm<sup>3</sup> *η<sup>b</sup>* Myofibroblast to fibroblast collagen production rate proportionality factor 2 -

*<sup>d</sup><sup>ρ</sup>* Collagen degradation rate per unit of cell density 7.59 <sup>×</sup> <sup>10</sup>−<sup>8</sup> cm3/cell day

Parameter Description Value Dimension *p*max Maximal cellular active stress per unit of ECM 10−<sup>4</sup> N g/cm2 cell *<sup>K</sup>*pas Volumetric stiffness moduli of the passive components of the cell <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> N g/cm2 cell *<sup>K</sup>*act Volumetric stiffness moduli of the active filaments of the cell 1.852 <sup>×</sup> <sup>10</sup>−<sup>5</sup> N g/cm2 cell

*τ<sup>d</sup>* Half-maximal mechanical enhancement fo fibroblast to myofibroblast differentation 10−<sup>5</sup> N g/cm2 cell *μ*<sup>1</sup> Undamaged skin shear viscosity 200 N day/cm2 *μ*<sup>2</sup> Undamaged skin bluk viscosity 200 N day/cm2 *E* Undamaged skin Young's modulus 33.4 N/cm2 *ν* Undamaged skin Poisson's ratio 0.3 *ξ* Myofibroblasts enhancement of traction per unit of fibroblasts density 10−<sup>3</sup> cm3/cell *<sup>R</sup><sup>τ</sup>* Traction inhibition collagen density <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> g/cm<sup>3</sup> *<sup>s</sup>* Dermis tethering factor <sup>5</sup> <sup>×</sup> 102 N/cm g

*θ*<sup>1</sup> Shortening strain of the contractile element −0.6 *θ*<sup>2</sup> Lengthening strain of the contractile element 0.5 -

Table 3. Parameters related to the mechanical behaviour of cells and the ECM.

*η<sup>d</sup>* Myofibroblast to fibroblast collagen degradation rate proportionality factor 2 -

*ζ* Myofibroblast to fibroblast chemical production rate proportionality factor 1 - <sup>Γ</sup> Half-maximal enhancement of net GF production <sup>×</sup>10−<sup>8</sup> g/cm<sup>3</sup> *dc* GF decay rate 0.693 day−<sup>1</sup>

Table 2. Parameters related to the collagen and growth factor equations.

ecm Initial GF concentration in the wound 10−<sup>8</sup> g/cm<sup>3</sup> *Dc* GF diffusion rate <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> cm2/day *kc* GF production rate per unit of cell density 7.5 <sup>×</sup> <sup>10</sup>−<sup>6</sup> cm3/cell day

*λ*0

*c*0


24 Will-be-set-by-IN-TECH

512 Tissue Regeneration – From Basic Biology to Clinical Application

Javierre, E., Moreo, P., Doblaré, M. and García-Aznar, J.M.: *Numerical modeling of a*

Vermolen, F.J. and Javierre, E.: *A suite of continuum models for different aspects in wound healing*,

Vermolen, F.J. and Javierre, E.: *Computer simulations from a finite-element model for wound*

Vermolen, F.J. and Javierre, E.: *A finite-element model for healing of cutaneous wounds combining*

Schugart, R.C., riedman, A., Zhao, R., Sen, C.K.: *Wound angiogenesis as a function of tissue oxygen*

Xue, C., Friedman, A., Sen, C.K.: *A mathematical model of ischemic cutaneous wounds*, Proceedings of the National Academy of Sciences USA, 106(39), 16783–16787, 2009. Graner, F., Glazier, J.: *Simulation of biological cell sorting using a two-dimensional extended Potts*

Merks, M.H., Koolwijk, P.: *Modeling morphogenesis in silico and in vitro: towards quantitative,*

Plank, M.J. and Sleeman, B.D.: *Lattice and non-lattice models of tumour angiogenesis*, Bull.

Vermolen, F.J. and Gefen, A.: *A semi-stochastic cell-based formalism to model the dynamics of*

Kumar, V., Abbas A. and Fausto, N.: *Robbins and Cotran Pathologic Basis of Disease*, Philadelphia:

Komori, M., Tomizawa, Y., Takada, K. and Ozaki, M.: *A single local application of recombinant*

*human basic fibroblast growth factor accelerates initial angiogenesis during wound healing*

engineering and biomaterials, New York: Springer-Verlag; 2009.

*contraction and closure*, Journal of Tissue Viability 2010; 19:43-53.

*contraction, angiogenesis and closure*, submitted, 2011.

*model*, Phys. Review Letters, *69*, 2013–2016, 1992.

Mathem. Biol., *66* (2004), 1785–1819.

*migration of cells in colonies*, to appear, 2011.

*in rabbit ear chamber*, Anesth Analg, *100* (2005), 830–834.

and Structures 2009; 46:3597-3606.

105(7), 2628–2633, 2008.

Elsevier Saunders; 2005.

149–171, 2009.

*mechano-chemical theory for wound contraction analysis*, International Journal of Solids

In: Bioengineering research of chronic wounds, studies in mechanobiology, tissue

*tension: a mathematical model*, Proceedings of the National Academy of Sciences USA,

*predictive, cell-based modeling*, Mathematical modeling of natural phenomena, *4(4)*,

### *Edited by Jamie Davies*

When most types of human tissue are damaged, they repair themselves by forming a scar - a mechanically strong 'patch' that restores structural integrity to the tissue without restoring physiological function. Much better, for a patient, would be likefor-like replacement of damaged tissue with something functionally equivalent: there is currently an intense international research effort focused on this goal. This timely book addresses key topics in tissue regeneration in a sequence of linked chapters, each written by world experts; understanding normal healing; sources of, and methods of using, stem cells; construction and use of scaffolds; and modelling and assessment of regeneration. The book is intended for an audience consisting of advanced students, and research and medical professionals.

Tissue Regeneration - From Basic Biology to Clinical Application

Tissue Regeneration

From Basic Biology to Clinical Application

*Edited by Jamie Davies*

Photo by pistolseven / iStock