Original Research ARTICLE
A Preliminary Computational Investigation Into the Flow of PEG in Rat Myocardial Tissue for Regenerative Therapy
- 1Department of Mechanical Engineering, University of Cape Town, Rondebosch, South Africa
- 2Wallenberg Research Centre, Stellenbosch Institute of Advanced Study, Stellenbosch University, Stellenbosch, South Africa
- 3Department of Mechanical Engineering, University College London, London, United Kingdom
- 4Cardiovascular Research Unit, Department of Surgery, University of Cape Town, Observatory, South Africa
- 5Division of Biomedical Engineering, Department of Human Biology, University of Cape Town, Observatory, South Africa
- 6Bioengineering Science Research Group, Engineering Sciences, Faculty of Engineering and the Environment, University of Southampton, Southampton, United Kingdom
Myocardial infarction (MI), a type of cardiovascular disease, affects a significant proportion of people around the world. Traditionally, non-communicable chronic diseases were largely associated with aging populations in higher income countries. It is now evident that low- to middle-income countries are also affected and in these settings, younger individuals are at high risk. Currently, interventions for MI prolong the time to heart failure. Regenerative medicine and stem cell therapy have the potential to mitigate the effects of MI and to significantly improve the quality of life for patients. The main drawback with these therapies is that many of the injected cells are lost due to the vigorous motion of the heart. Great effort has been directed toward the development of scaffolds which can be injected alongside stem cells, in an attempt to improve retention and cell engraftment. In some cases, the scaffold alone has been seen to improve heart function. This study focuses on a synthetic polyethylene glycol (PEG) based hydrogel which is injected into the heart to improve left ventricular function following MI. Many studies in literature characterize PEG as a Newtonian fluid within a specified shear rate range, on the macroscale. The aim of the study is to characterize the flow of a 20 kDa PEG on the microscale, where the behavior is likely to deviate from macroscale flow patterns. Micro particle image velocimetry (μPIV) is used to observe flow behavior in microchannels, representing the gaps in myocardial tissue. The fluid exhibits non-Newtonian, shear-thinning behavior at this scale. Idealized two-dimensional computational fluid dynamics (CFD) models of PEG flow in microchannels are then developed and validated using the μPIV study. The validated computational model is applied to a realistic, microscopy-derived myocardial tissue model. From the realistic tissue reconstruction, it is evident that the myocardial flow region plays an important role in the distribution of PEG, and therefore, in the retention of material.
The prevalence of non-communicable diseases has risen significantly in the last few decades, and this has affected societies in a myriad of ways. In low- and middle-income countries, cardiovascular and other non-communicable diseases tend to affect younger people who are economically active, and this has a negative impact on the general well-being of these nations (1, 2). Recognition of these changing trends has led to various shifts in global health policy (3, 4). From a healthcare management perspective, greater focus has been placed on prevention of lifestyle diseases, but also on improved understanding of disease progression. This has created an impetus for the development of appropriate solutions for long term management or mitigation of disease. In the last decade, regenerative medicine has gained a lot of attention as a potential long-term solution for a broad range of diseases. Included in these conditions of interest is myocardial infarction (MI), more commonly known as heart attack (5).
The end result of MI is the loss of a large number of cardiomyocytes and deposition of scar tissue at the site of injury. This leads to a decrease in left ventricular function, placing strain on the rest of the heart to keep supplying the body with adequate volumes of blood at the correct rate. Existing treatments for myocardial infarction focus on limiting the amount of damage that the tissue suffers but ultimately, they just serve to prolong the time to heart failure (6, 7). While it is highly unlikely that the injured myocardium can ever return to its original state, some functional benefit has been gained from injecting stem cells into the injured part of the heart (8). The exact mechanism of healing following stem cell engraftment is still unclear. Some studies have shown that the paracrine effects of stem cell associated growth factors contribute significantly to healing (9–12).
The main challenge with the efficacy of stem cell treatment in myocardial infarction, in particular, is the low retention of cells following injection (13). Owing to the vigorous action of the heart, a number of injected cells are transported via the systemic circulation to other organs of the body, while others are ejected via the needle track (14–16). Factors such as delivery route and injection timing have been investigated as possible contributors to retention (17–22). Even with the massive losses of injected cells, improved left ventricular function has been observed. This has encouraged development of scaffolds and biomaterials which can be injected with the stem cells, aiding the engraftment process. In some instances, the biomaterials have been seen to improve ventricular function without the stem cells, leading to numerous studies attempting to optimize their properties (23). Hydrogels are a type of scaffold that has been widely explored for these purposes and can be derived from natural or synthetic materials.
Polyethylene glycol (PEG) gel and its modified forms have been used as an injectate for myocardial infarction therapy development (24–28). Synthetic PEG gels can be readily customized and optimized for specific purposes. When a suitably functionalised PEG is prepared in solution and combined with a cross-linker, it can be injected into the heart in a fluid state and then spontaneously gelled over time, in situ. A key observation that has been made about the gel following injection is that it distributes itself differently in the myocardial tissue, depending on the injection timing (24). Less gel is retained if injection happens immediately after infarction while delayed delivery results in improved retention and functional benefit. Such observations have led to a greater desire to understand the mechanism of PEG distribution and retention once injected.
To this end, a number of computational models have been developed in an effort to understand the behavior and effect of injectable scaffolds (29–37). Most of these models have focused on biomaterials in their gelled state and have examined how the presence of the gel in the myocardial tissue affects stress response in the wall. Given that the greatest losses are likely to occur before the gelled state is reached, interest has been expressed in the behavior of the gel while still in a liquid state. While there are currently no in silico studies examining the behavior of liquid PEG in myocardial tissue, there are a vast number of experimental studies that have examined the behavior of PEG in a wide variety of circumstances. These often consider the modified forms of PEG or the combined response of the gel and stem cells (38–41). The main challenge is that only a very small subset of these studies focus on PEG in its pre-altered form, hence making it difficult to state its fundamental rheological behavior. A number of studies focusing on the PEG monomer have examined how variables such as pH, molecular weight, and temperature alter rheological behavior (42–48).
For most of the studies focusing solely on PEG, the material has been shown to have constant viscosity within a specified shear rate range, and this value varies according to the molecular weight of the gel and the testing conditions. Dontula et al. used PEG 8000 as a reference in a study examining the behavior of PEG solutions (42). The reference solution exhibited Newtonian behavior up to a shear rate of 100 s−1. In a study that examined the effect of thermal history on the behavior of PEG (3, 6, and 10 kDa), Baldursdottir et al. found that viscosity decreased dramatically within the first 1°C of temperature change and then stabilized at a constant value for each different molecular weight (44). The experiments were all conducted at the same shear rate (0.1 s−1), hence it was not possible to conclude how the viscosity would change with shear rate. It was evident, however, that following the initial drop in temperature, viscosity generally remained constant. Azri et al. characterized the linear viscoelastic behavior of low molecular weight PEG by studying the rheological response of 0.6 kDa PEG over a frequency range of 0.16–16 Hz (46). The experiments were conducted at different temperatures and it was found that all samples exhibited Newtonian behavior over the entire range of frequencies. In a study using 4 kDa PEG, Brikov et al. also found that viscosity is constant until 1,000 s−1 (47). Contrary to this general trend, Liu et al. found that PEG was shear-thinning over a shear rate range of 2,000 s−1 (48). In their study, they considered the rheological character of PEG over a range of molecular weights (0.4, 0.6, 1, 6, 10, and 20 kDa). For all the samples, viscosity was seen to decrease as shear rate increased.
The wide range of testing conditions and methods employed for characterizing the rheological behavior of PEG make it difficult to draw exact conclusions. For example, some of the studies make use of a viscometer, which produces information about viscosity, but not about the overall flow behavior of a fluid (43, 45, 49). In addition, available literature on PEG tends to focus on a particular molecular weight or temperature, which is of interest for further studies. The assumption that the relative simplicity of the material would have resulted in complete characterization does not seem hold, and the lack of standard testing conditions (e.g., concentrations, temperatures, and pH) makes it difficult to compare results. Furthermore, most of the rheological techniques in these studies characterize fluid behavior on the macroscale. This may result in incorrect characterization of dynamics that arise on the microscale (50–52). This is particularly true in the case of polymers, which experience restrictions on conformation in near-wall regions and behave slightly differently in microflows (53, 54).
In this study, we aim to characterize and study the flow behavior of a 20 kDa PEG solution in myocardial tissue, after injection but prior to gelation, when the greatest losses are likely to occur. The vast majority of results presented from rheological tests, in literature, state that PEG behaves as a Newtonian fluid within a specific shear rate range. Its microscale flow behavior, where viscous forces are likely to be dominant, is unknown and is an important factor to consider when modeling flows in myocardial tissue. Micro particle image velocimetry (μPIV) was used to observe the flow profiles of PEG solution in a microchannel, which represents the gaps in myocardial tissue, and is comparable in size. Even though extensive remodeling takes place within the myocardial tissue following infarction, the early stage post-infarct, which is considered in this model, is adequately represented by this microchannel. The data obtained from these μPIV experiments was used to characterize the behavior of the fluid on the microscale. An idealized two-dimensional computational fluid dynamics (CFD) model, replicating the experimental microchannel, was then developed and validated using the μPIV data. Even though the geometry in question is fairly simple, validation is important as there are no analytical solutions for the flow of non-Newtonian fluids in non-circular geometries. The idealized CFD model was then extended to a three-dimensional model of realistic myocardial tissue derived from confocal microscopy. The myocardial tissue resembles normal tissue, and, like the microchannel, reflects the early stage post-infarct. As a proof of concept, this model was used to investigate whether or not the myocardial flow region had a significant impact on gel transport.
Materials and Methods
The study comprised two main parts. The in vitro study sought to establish the flow behavior of PEG solution in microchannels, using μPIV. The in silico study, which used some of the data from the in vitro study, then examined flow in idealized and realistic geometries.
In vitro Study
In order to characterize the flow of PEG in microchannels representing the gaps in myocardial tissue, a solution of PEG had to be prepared. This was then perfused through a microchannel, and the flow was imaged and quantified using a μPIV setup. Finally, a set of equations was used to characterize flow behavior and extract key parameters. Details of each stage are given below.
Polyethylene Glycol Solution Preparation
PEG solution was prepared using a protocol adapted from the work of Kadner et al. (24). A 20 kDa, 8-arm hydroxyl-terminated PEG (20PEG-8OH, Nektar) was used for the study. A 10% m/v solution of PEG was prepared by dissolving 1 g of powdered PEG-OH in 10 ml PBS. Unlike other injectate studies where a cross-linker is added for gelation, the solution was left as is to enable observations in the liquid state, when the greatest losses are most likely to occur in vivo.
The flow of the solution was observed in a polydimethylsiloxane (PDMS) microchannel with a square cross-section of 50 × 50 μm and length of 15 mm. The dimensions of the microchannel approximated the gaps through which the solution would flow in myocardial tissue. The flow rate of the solutions was regulated using a custom made pneumatic microfluidic flow controller as described in Sherwood et al. (55). The microchannel was placed on an inverted microscope (DMILM: Leica). The flow was seeded with 1 μm density-neutral fluorescent microspheres (Nilon Red, Invitrogen) and illuminated by means of a pulsed Nd:YAG laser (wavelength 532 nm, New Wave, USA). The focal plane was position at the central plane of the channel and 60 pairs of images were acquired by a CCD camera (Hamamatsu, UK:C8484-05C) with 10x magnification and sampling rate of 6 Hz. The time interval between the pairs of images was 0.5 ms. Image acquisition was controlled via LabVIEW (National Instruments) and images were pre-processed with MATLAB. A schematic of the μPIV system is illustrated in Figure 1. The velocity was obtained from the acquired images using multipass PIV cross-correlation algorithms implemented in open source software JPIV. Three pass ensemble average correlation was utilized with interrogation windows of 8 pixels height and widths of 32, 16, and 16 pixels, respectively, with 50% overlap resulting in a vector field with a spatial resolution of 2.4 μm across the flow. A region of interest (ROI) three channel widths long and located at the middle of the microchannel, i.e., 7.5 mm downstream of the inlet, was selected to estimate the velocity profiles across the channel. The distance from the wall to the first measurement point was 2 μm. The nominal shear rate, defined as the ratio of the average velocity to channel width, ranged from 8 to 132 s−1 in this study. The average velocity (and the flowrate) was determined by integrating the measured velocity profiles. The estimated flow rates (ranging from 2.19 × 10−8 to 3.42 × 10−7 m3s−1) were used as input in the CFD simulations for the simplified channel geometry. Thus, any effects due to depth of correlation (equal to 36 μm in this set up) will be minimized or cancel out when making experimental and CFD comparisons.
Characterizing Fluid Behavior From μPIV Velocity Profiles
The data obtained from the μPIV study was analyzed and fitted using equations described in the paragraphs that follow. From this data, the rheological behavior of the fluid was deduced. Details of the steps are given in the paragraphs that follow.
The measured velocities across the microchannel were fitted to an experimental data curve derived by Sherwood et al. using MATLAB's non-linear fitting function (55). This experimental curve corrects for irregularities and asymmetries which are present in the measured data. The curve used for fitting is given by Equation (1),
where vexperimental is velocity for the experimental curve, v is a weighting parameter related to average velocity, m is an exponent that describes the bluntness of the velocity profile, is the y co-ordinate normalized by width (y/w) and ranges from −0.5 to 0.5, and v0 is the wall slip velocity. Wall slip was included in the fitting as PIV readings could not be accurately obtained near the wall; hence the first velocity measurement was taken 4 pixels from the wall. Furthermore, it is important to consider wall slip in polymer solutions, which have features that deviate from the no-slip condition. The slip layer, which is influenced by polymer concentration and particle size, can be >2 μm (53). Often, the Equation (1) is a simplification of the analytical velocity profile for laminar Newtonian flow in rectangular channels (56). The two-dimensional analytical velocity profile is given in Equation (2),
where vanalytical is the analytical velocity for a given y position, h is the height of the channel, Δp is the pressure gradient across the channel length, μ is the dynamic viscosity, L is the length of the channel, and w is the channel width. The deviation of the measured velocity profiles from the analytical solution provides a measure of the extent of shear thinning of the fluid as the latter tend to produce blunter profiles.
Once the experimental curves were obtained, they were normalized with the average velocity so that the power law index for each case could be retrieved. It is important to note that there is no analytical solution for the flow of a power law fluid in a rectangular duct (52). Fortunately, the power law index is independent of the cross-sectional shape of the duct (57). For this reason, the power law index can be calculated using the analytical solution for a power law fluid in a circular duct. The power law index for each case was then determined using Equation (3),
where vnormalized is the experimental velocity normalized by average velocity, n is the power law index, vaverageN is the average velocity calculated for the fitting and v0N is the slip velocity calculated for the fitting. MATLAB's curve fitting tool was used to calculate n and vaverageN, and the Levenberg-Marquardt algorithm was selected in the non-linear least squares method. The value of n indicates whether the fluid is Newtonian or non-Newtonian, with n = 1 for a Newtonian fluid, n > 1 for a shear-thickening fluid, and n < 1 for a shear-thinning fluid.
Differential capillary viscometry, a technique which has been used to determine viscosity in microchannels, presents a method for deducing viscosity from flow, and a more detailed description can be found in Wunderlich and Bausch (52). By integrating the power law index with respect to shear rate for each case, the viscosity-shear rate curve can be recovered. The viscosity for each case is calculated using Equation (4),
where is the shear rate. Once the viscosities are known, the power law relation is then used to calculate the consistency index, K, using a log-log relation for viscosity and shear rate. This fitting also gives an ideal power law index, which is used for the realistic cases. The power law equation is given by Equation (5),
where K is the consistency index. Even though polymers tend to be better characterized by other fitting models, a simplified version of the power law equation is used to estimate parameters in an attempt to minimize error introduced by the estimation of too many parameters.
In silico Study
The CFD study made use of two different geometries to validate and observe flow behavior in myocardial tissue. An idealized channel was used to compare computational results with experimental data. Confocal microscopy images of rat myocardium were reconstructed to develop a model based on realistic representation of myocardial tissue. This was used to observe flow in a more complex setup and to investigate the impact of myocardial flow region.
An idealized channel was modeled and implemented in two-dimensions. The flow domain has a width of 50 μm and a length of 15 mm, as illustrated in Figure 2. These dimensions were chosen to match the dimensions of the microchannel used in the experimental study, and to enable validation of the computational study. Geometries were developed using DesignModeler and meshes were generated in ANSYS Meshing (ANSYS Inc., Lebanon, NH).
A grid independence study was conducted to determine the appropriate resolution level for the idealized simulations in this study. Grids of increasing refinement were defined, with element sizes of 6.25 × 10−6, 4.6875 × 10−6, 3.7594 × 10−6, and 3.125 × 10−6 m. The percentage error between the different meshes was found to be between 1 and 2%, and all results were computed on a mesh of element size 3.125 × 10−6 m. All the cells in the computational domain were the same size.
Realistic Myocardial Tissue
The realistic, three-dimensional geometry was obtained from confocal microscopy images of a rat's heart, illustrated in Figure 3 (58). Confocal microscopy is an optical imaging technique used for increasing optical resolution and contrast of a micrograph. It captures two-dimensional images at different depths of a sample and allows the reconstruction of these images into three dimensional structures. This geometry was used to investigate the manner in which three different myocardial flow regions, labeled 1, 2, and 3 in Figure 3, affect PEG dispersion in the tissue. The length of the tissue from endocardium to epicardium is ~4.25 mm.
Figure 3. Confocal microscopy image of left ventricular tissue slice from a rat. Adjusted with permission from Sands et al. (58).
Figure 4 illustrates a magnified view of heart tissue. Various regions of interest can be observed in the image. A slice taken from the mid-wall region of a rat's left ventricle is shown in Figure 4A. This image is captured under a 25x magnification lens. The region labeled (P) indicates a cleavage plane, which is the interstitial space located between and around cardiomyocyte bundles, while (V) illustrates the microvasculature structure. A 63x magnification of the same tissue is illustrated in Figure 4B. The myocytes are indicated by (M) and the collagen chords by (C). In the computational model developed, it is assumed that flow of the hydrogel takes places in the cleavage planes found in the tissue.
Figure 4. (A) Slice taken from a mid-wall region of a rat's left ventricle analyzed under a 25x magnification lens. (P) indicates the cleavage planes and (V) the vasculature structure. (B) 63x magnification of the same tissue in (A). Here, the myocytes can be seen (M) and collagen chords in (C). Picrosirius red was used for the staining of collagen. Adjusted with permission from Young et al. (59).
Image reconstruction was carried out in Simpleware ScanIP (Simpleware, Exeter, UK). Stacks of high resolution and low resolution images were obtained from confocal microscopy. Using a thresholding technique, distinction was made between cleavage planes and solid myocardium, enabling reconstruction of the tissue block. Manual segmentation was employed and in some instances, myocardial material was removed to enable continuity between the cleavage planes. This was done to ensure that the fluid flow equations would not be violated. To examine the flow through the myocardial tissue, blocks of dimensions 660 × 560 × 440 μm were chosen as the region of interest, as illustrated in Figure 3. The blocks were taken as close to the epicardium as possible, where injection is likely to take place. Tetrahedral meshes of varying refinement were created using the Simpleware FE module. Figure 5 illustrates a sample mesh that was created. The gray zone represents the cleavage plane and is defined as a fluid region through which PEG can flow. The white region represents the tissue and is modeled as a rigid wall.
Figure 5. Tetrahedral mesh of heart tissue created in Simpleware. The gray zones represent areas where fluid can flow while white zones represent solid tissue.
A grid independence study was carried out for the realistic geometry. In this instance, it was slightly harder to control element sizes as tetrahedral elements were used. The coarsest mesh had 274 477 cells, followed by a medium mesh of 425 136 cells and a fine mesh of 730 799 cells. The error between the medium and fine mesh was 1.5% and results were computed on the medium mesh.
In silico Settings and Parameters
For all the above-mentioned cases, the Navier-Stokes equations were used to describe the flow of the fluid. The continuity equation is described by Equation (6) and the momentum equation by Equation (7),
where ρ is density, v is the velocity vector, t is time, and P is pressure.
The computational frameworks were implemented in ANSYS Workbench 19.2 (ANSYS Inc., Lebanon, NH). Flow simulations were computed in FLUENT and CFD-Post was used for post-processing. PEG is modeled as a non-Newtonian, incompressible fluid with density ρ = 1,200 kg.m−3 (60). For viscosity, a non-Newtonian power law is used and is described in Equation (5). For this study, the gel is modeled in its most fluid state, prior to the start of gelation. In reality, PEG gel is a viscoelastic material with a non-Newtonian viscosity that changes as the solution gels. All the simulations were run under steady state conditions in order to match the conditions employed for obtaining the μPIV data. A mass flow rate was specified at the inlet boundary and a zero pressure condition at the outlet boundary. For the idealized geometries, slip was specified at the walls as the fluid in question is a polymer solution flowing in a small channel (53, 54). The slip for the realistic geometry could not be determined hence a no-slip condition was employed for those cases. With the exception of density, all other parameters were calculated from the experimental study and are presented in the results section. The viscosity model values and the mass flow inlet for the realistic tissue model are also calculated from experimental values and are presented in the results section. For most of the μPIV results, errors are minimized by normalization. For the purposes of capturing the varying flow rate in the CFD study, it was necessary to use the absolute μPIV variables.
SIMPLE algorithm was employed for pressure-velocity coupling. A least squares cell based scheme was used for the gradient, and both pressure and momentum were solved with a second-order scheme. An algebraic multigrid solver was used for calculations. In order to validate the computational results, comparisons were made with those obtained from the μPIV experiments.
In this section, the results of the μPIV and computational studies are presented. The comparison between experimental and computational results is presented in the computational results section.
In vitro Study
Figure 6 illustrates the experimental curves following fitting of the original velocity data. The mass flow rate increases from case 01 to case 13, with case 01 reaching a maximum velocity of ~5.8 × 10−4 m.s−1 and case 13 exhibiting a maximum velocity of 9.2 × 10−3 m.s−1. The average velocity for case 01 is 4.4 × 10−4 m.s−1 and 6.9 × 10−3 m.s−1 for case 13. The ratio of the maximum to average velocity is ~0.75 for all the cases shown in the figure. The wall slip for each case, which also happens to be the minimum velocity, follows a similar trend to the other variables, with a value of 1.5 × 10−4 m.s−1 for case 01 and 2.2 × 10−3 m.s−1 for case 13.
Figure 6. Experimentally measured flow velocities across the channel width for a range of flow rates (μPIV study). EXP01 is the lowest flow rate and EXP13 the highest.
The experimental curves were normalized using the appropriate average velocity for each case, as illustrated in Figure 7. The area under each curve is 1, the maximum value is ~1.3 and the minimum value is ~0.3. The solid black line represents the analytical Newtonian profile, and it is evident that the experimental curves are blunter than the analytical solution implying shear thinning behavior.
Figure 7. Normalized, experimentally derived, flow velocities across channel width for a range of flow rates. The solid black line represents the analytical profile for a Newtonian fluid.
Figure 8 illustrates the fitting of the normalized data to the power law curve. The example chosen for this specific image is Case 13 and the power law index, n, was 0.9389, indicating slight shear-thinning behavior. Other parameters were also calculated for the fitting, with vaverageN = 0.5198 and v0N = 0.3186. Fitting was carried out for each case.
The calculated values obtained from the experimental data and the various fittings are detailed in Table 1. The shear rate ranges from 8 s−1 for Case 01 to 132 s−1 for Case 13. The flow rate ranges from 2.19 × 10−8 to 3.42 × 10−7 m3s−1. The power law index for most of the cases indicates shear-thinning behavior, and ranges from 0.793 to 0.943. Even with the anomaly presented by case 10, where the value indicates Newtonian/shear-thickening behavior, the fluid is assumed to be shear-thinning. Excluding case 10, viscosity varies from 0.434 to 0.887 kg.m−1s−1. In theory, the power law index remains constant for a given fluid. The power law indices for each of the flow rates are determined by fitting the normalized velocity data to the power law curve for a circular duct, on a case by case basis, as shown in Figure 8. Some variations are observed from case to case, and this is reflected by the slight variation in bluntness for each of the normalized velocity profiles, as illustrated in Figure 7.
Finally, the consistency index for the fluid was found to be 1.0 kg.sn−2.m−1 and the overall power law index of the fluid is 0.883. These two parameters were used for the viscosity model of the realistic tissue model, and the flow rate for case 13 was used as the inlet mass flow rate for the same model.
In silico Study
Figure 9 illustrates a comparison between the experimental and the CFD curves for the same parameters. The fluid property results presented in the in vitro study were applied to the computational model. The solid lines represent the CFD data, the dots represent the experimental values and the dashed lines represent the 99% confidence interval for those experimental values. Most of the CFD curves fit within the 99% confidence interval of the experimental values. Case 8 and 11 stand out as they have profiles that are significantly blunter than the experimental values, resulting in some of the values falling outside the 99% confidence interval. This discrepancy decreases with decreasing flow rate. Case 01, which has the lowest flow rate, has a maximum experimental velocity of 5.8 × 10−4 m.s−1 and a maximum CFD velocity of 5.7 × 10−4 m.s−1. The mean experimental and CFD values are 4.4 × 10−4 and 4.3 × 10−4 m.s−1, respectively. Case 13, which has the highest flow rate, has a maximum velocity of 9.2 × 10−3 m.s−1 for the experimental case, while the maximum CFD velocity is 9.0 × 10−3 m.s−1. The mean experimental value for the same case is 6.9 × 10−3 m.s−1 and the mean CFD value is 6.7 × 10−3 m.s−1.
Figure 9. A comparison between the experimentally measured and computationally predicted flow velocities across the channel width at the mid-plane, for a range of flow rates. The solid lines represent CFD values, the dots represent experimental values and the dashed lines indicate the 99% confidence interval for those experimental values.
A plot of the normalized values for the CFD and experimental curves is illustrated in Figure 10. The curves all follow a similar trend, and all the values fall within the normalized 99% confidence interval of Case 13. The area under each curve is 1. The maximum normalized velocity values for the CFD plots lie between 1.31 and 1.37. The minimum normalized velocity values, which correspond to the slip velocity, range from 0.27 to 0.35, for the CFD plots. All the profiles are blunter than the Newtonian profile, indicating shear-thinning behavior.
Figure 10. Comparison of normalized experimental and two-dimensional CFD data. The solid lines represent CFD values, the dots represent experimental values, the dashed lines indicate the 99% confidence interval for those experimental values and the black solid line represents the analytical, Newtonian solution.
The non-Newtonian flow properties for the overall fluid were applied to different myocardial flow regions, as illustrated in Figure 11. It is clear, from the illustration, that the myocardial flow region has an impact on the manner in which the fluid passes through the tissue. Sections with slightly wider cleavage planes, such as Block 2 and Block 3, tend to achieve lower maximum flow velocities. The maximum velocity achieved for Block 2 is 2.15 × 10−3 m.s−1, and this is recorded on the outlet plane. Block 3, which is also has relatively wide cleavage planes achieves a maximum velocity of 4.55 × 10−3 m.s−1, however this value is not visualized on any of the planes illustrated in the diagram. The maximum value for block 1 is much higher than those measured on the selected planes. As such, the inlet, mid-plane, and outlet appear to have very little activity. The maximum flow velocity for Block 1 is 1.21 × 10−1 m.s−1 and this value is not visualized on any of the selected planes. This value is approximately two orders of magnitude higher than the maximum values for Blocks 2 and 3.
Figure 11. Velocity contour plots for each myocardial flow region, labeled Block 1, Block 2, and Block 3. The block names correspond to the sections illustrated in Figure 4. Velocity contours are plotted at the inlet, mid-plane, and outlet of each block. The red line indicates the centerline through the XY mid-plane, for which the velocity, shear rate, and viscosity are illustrated in Figure 12.
The velocity profiles along the xy mid-plane of each block, indicated by the red lines in Figure 11, are illustrated in Figure 12. The values reflected here are indicative of the flow velocities, shear rates, and viscosities experienced on the mid-plane, and may not necessarily reflect the flows experienced in the overall block. The same mass flow rate is specified at the inlet of each block. Even though the contour plots of Block 1 show very little activity in Figures 11, 12 clearly illustrates that the velocity values experienced on the mid-plane are fairly similar to those experienced on the mid-plane of Block 2. Block 1 has a maximum flow velocity of 5.60 × 10−4 m.s−1, while Block 2 has a peak flow velocity of 5.20 × 10−4 m.s−1. A slightly lower peak is also present for Block 2. Of the three blocks, Block 3 exhibits the highest peak, which has a value of 1.10 × 10−4 m.s−1.
Figure 12. Velocity profiles for the different blocks taken from a line on the xy mid-plane, as indicated in Figure 11.
Figure 13, which illustrates the shear rate within the channels, clearly demonstrates the range of conditions that the fluid would be exposed to within a single channel. The non-Newtonian behavior of the fluid becomes particularly important under these circumstances, as variation in shear rate would result in corresponding variations in viscosity. For a shear thinning-fluid, an increase in shear rate would result in a decrease in viscosity, thus decreasing the fluid's resistance to flow. For block 1, the fluid experiences a shear rate range of 1–36 s−1, within a single channel. The ranges for blocks 2 and 3 are 0.7–28 s−1 and 2–38 s−1, respectively. These results reflect the heterogeneity along a single line, on a single plane. Calculation of shear rates across the entire block would yield even greater variation.
Discussion and Conclusions
Existing computational models of hydrogel injectates in myocardial tissue have examined injectates in the gelled state and have elucidated the impact of these gels on stresses in the myocardium (29–37). These studies have been useful in elucidating some of the mechanisms that contribute to the long-term efficacy of hydrogels. One of the main gaps in the literature is retention studies examining the manner in which fluids flow through myocardial tissue, prior to gelation. The work presented in this study makes a contribution toward understanding how fluids flow in myocardial tissue, and thereby how fluids are lost from the myocardium, and presents a framework for further studies of this nature. This complements observations such as those made by Van den Akker et al. who were able to visualize, in real-time, the wash out of stem cells via the venous circulation (61).
Dontula et al. who used 8 kDa PEG, found that PEG behaved as a Newtonian fluid up to 100 s−1 (42). Beyond this value, the viscosity of the solution was seen to increase, indicating shear-thickening behavior. Brikov et al. who used 4 kDa PEG, observed a similar phenomenon, although they found that the viscosity began increasing at a higher shear rate (47). These studies used standard, macroscopic, rheological techniques to reach their conclusions. In this study, there was interest in characterizing the gel on the microscale, where its behavior could deviate from macroscale observations. This is particularly true in the case of polymers, which experience restrictions on conformation in near-wall regions and behave slightly differently in microflows (53, 54). μPIV and CFD were used to examine the flow of PEG solution on the microscale, in channels resembling gaps found in heart tissue. A method presented in this paper, where the rheological parameters were derived from μPIV flow results, revealed that PEG solution behaves as a shear-thinning fluid on the microscale for shear rates between 8 and 129 s−1. These findings are in agreement with the work of Liu et al. who found that PEG was shear-thinning up until a shear rate of 2,000 s−1 (48). This was found to be true for PEG solutions of molecular weights ranging from 0.4 to 20 kDa. It is interesting to note that even lower molecular weight PEG solutions/gels exhibited different properties on the nanoscale to what was previously reported. There are various factors which may give rise to this Newtonian/non-Newtonian discrepancy seen in literature, including the testing conditions, molecular weight of PEG gel and the concentrations of PEG in the solutions. A more thorough investigation spanning spatial scales, testing methods, PEG molecular weights, and concentrations would be beneficial in determining how the material behaves for a range of different conditions.
The experimental results were useful not only for deducing the rheological properties of the fluid on the microscale, but also for providing data for input into the computational model. First, an idealized, two-dimensional CFD model was developed to enable meaningful comparison with experimental data, and for validation purposes. There was good agreement between absolute and normalized values, indicating that the methodology for deducing the rheological properties was effective. The computational framework was then implemented in realistic, confocal microscopy-derived microstructural geometry of rat myocardial tissue. This model was used to examine the effect of myocardial flow region on flow patterns within the tissue. It was evident, particularly from the shear rate results, that the fluid is exposed to a range of shear rates within individual channels, and this would lead to variation in viscosity values. Even though the channel widths were comparable for the different blocks, the distribution of shear rates within each channel is unique, hence the myocardial flow region does have an impact on fluid distribution. A number of in vitro studies have examined how different delivery routes affect retention (23). It was clear that the exact site of injection had an impact on fluid distribution. The channels in the sections explored in this study were of a comparable size, but even with that, differences in flow speed could be seen. Future studies could explore the flow of PEG gels in more complex microfluidic geometries with varying channel sizes mimicking the tissue microstructure (62). An abundance of narrow cleavage planes is likely to result in higher flow velocities and therefore, greater gel losses. Wider cleavage planes are likely to result in slow flow velocities, thereby encouraging more retention. This hypothesis could support the observations in a study by Kadner et al. where narrow channels characteristic of healthy tissue (and injection immediately after infarction) are likely to result in greater gel loss (24).
At this stage, there is no clarity regarding the exact shear rates that a hydrogel would be exposed to in vivo. It is highly likely that the shear rates will vary greatly over the cardiac cycle, and that some of values will be quite high. If PEG is indeed shear-thinning over a wide range of shear rates, the fluid would become less viscous as the shear rate increased, resulting in greater flow velocities and increased venous drainage. While this cannot be confirmed with the current model, the results presented here give clues as to why the retention of hydrogels has proven challenging to address. When the hydrogel is injected into the myocardium in a fairly liquid state, the solution is exposed to extensional flow as it moves through the narrow needle. Once the external pressure is removed and the gel moves through the myocardial tissue, the environment becomes far less predictable and the behavior of the solution prior to gelation is then harder to understand. The complex myocardial environment is likely to expose the injectate to a wide range of shear rates and may result in changes in the fluid's behavior, particularly on the microscale, where viscous forces tend to be dominant. The characterization of PEG solution over a wider range of shear rates, on the microscale, would therefore prove beneficial. In addition, a more realistic wall model which can account for the elastic effects of myocardial tissue, along with oscillatory shear experiments of the gel, would yield a more accurate representation of flow. Finally, a model which can account for the gelation process would enable verification of the period when the greatest losses occur.
The datasets generated for this study are available on request to the corresponding author.
MN contributed to the design of the work, acquisition, analysis, interpretation of data, drafting and critical revision of the manuscript. AP and SB contributed to the acquisition, analysis and interpretation of data and critical revision of the manuscript. JK, AL, JM, and LS contributed to the acquisition of data and to critical revision of the manuscript. DB contributed to the design of the work and to critical revision of the manuscript. ND and TF contributed to the conception and design of the work and to critical revision of the manuscript. All authors provide approval for publication and agree to be accountable for all aspects of the work.
This research was supported financially by the National Research Foundation (NRF grant UID92531 to TF and SFP13080926816 to MN) of South Africa. Any opinion, findings and conclusions or recommendations expressed in this publication are those of the authors and therefore the NRF does not accept any liability in regard thereto. AP was supported by the UCL Doctoral Training Programme in Medical Device Innovation. This programme was funded by University College London (UCL), EPSRC, National Institute for Health Research (NIHR) Biomedical Research Centres at University College London Hospitals, Great Ormond Street Hospital, Moorfields Eye Hospital NHS Foundation Trust and UCL Ophthalmology.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Mark Trew, who is based at the Cardiac Structure and Imaging Group at the Auckland Bioengineering Institute, University of Auckland, New Zealand, for providing confocal imaging data of rat myocardium.
5. Zhu Y, Matsumura Y, Wagner WR. Ventricular wall biomaterial injection therapy after myocardial infarction: advances in material design, mechanistic insight and early clinical experiences. Biomaterials. (2017) 129:37–53. doi: 10.1016/j.biomaterials.2017.02.032
6. Ravichandran R, Venugopal JR, Sundarrajan S, Mukherjee S, Sridhar R, Ramakrishna S. Minimally invasive injectable short nanofibers of poly (glycerol sebacate) for cardiac tissue engineering. Nanotechnology. (2012) 23:385102. doi: 10.1088/0957-4484/23/38/385102
9. Jiang Z-S, Srisakuldee W, Soulet F, Bouche G, Kardami E. Non-angiogenic FGF-2 protects the ischemic heart from injury, in the presence or absence of reperfusion. Cardiovasc Res. (2004) 62:154–66. doi: 10.1016/j.cardiores.2004.01.009
10. Segers VFM, Tokunou T, Higgins LJ, Macgillivray C, Lee RT, Macgillivray C, et al. Local delivery of protease-resistant stromal cell derived factor-1 for stem cell recruitment after myocardial infarction. Circulation. (2007) 116:1683–92. doi: 10.1161/CIRCULATIONAHA.107.718718
11. Laguens R, Meckert PC, Janavel GV, Del Valle H, Lascano E, Negroni J, et al. Entrance in mitosis of adult cardiomyocytes in ischemic pig hearts after plasmid-mediated rhVEGF165 gene transfer. Gene Ther. (2002) 9:1676–81. doi: 10.1038/sj.gt.3301844
12. Dries JL, Kent SD, Virag JAI. Intramyocardial administration of chimeric ephrinA1-Fc promotes tissue salvage following myocardial infarction in mice. J Physiol. (2011) 589:1725–40. doi: 10.1113/jphysiol.2010.202366
15. Li S-H, Lai TYY, Sun Z, Han M, Moriyama E, Wilson B, et al. Tracking cardiac engraftment and distribution of implanted bone marrow cells : comparing intra-aortic, intravenous, and intramyocardial delivery. J Thorac Cardiovasc Surg. (2009) 137:1225–1233.e1. doi: 10.1016/j.jtcvs.2008.11.001
16. Lyngbaek S, Ripa RS, Haack-Sorensen M, Cortsen A, Kragh L, Andersen CB, et al. Serial in vivo imaging of the porcine heart after percutaneous, intramyocardially injected 111In-labeled human mesenchymal stromal cells. Int J Cardiovasc Imaging. (2010) 26:273–84. doi: 10.1007/s10554-009-9532-4
17. Hou D, Youssef EA, Brinton TJ, Zhang P, Rogers P, Price ET, et al. Radiolabeled cell distribution after intramyocardial, intracoronary, and interstitial retrograde coronary. Circulation. (2005) 112:150–6. doi: 10.1161/CIRCULATIONAHA.104.526749
18. George JC, Goldberg J, Joseph M, Abdulhameed N, Crist J, Das H, et al. Transvenous intramyocardial cellular delivery increases retention in comparison to intracoronary delivery in a porcine model of acute myocardial infarction. Cell Ther. (2008) 21:424–31. doi: 10.1111/j.1540-8183.2008.00390.x
19. Hale SL, Dai W, Dow JS, Kloner RA. Mesenchymal stem cell administration at coronary artery reperfusion in the rat by two delivery routes: a quantitative assessment. Life Sci. (2008) 83:511–5. doi: 10.1016/j.lfs.2008.07.020
20. Mitchell AJ, Sabondjian E, Sykes J, Deans L, Zhu W, Lu X, et al. Comparison of initial cell retention and clearance kinetics after subendocardial or subepicardial injections of endothelial progenitor cells in a canine myocardial infarction model. J Nucl Med. (2010) 51:413–7. doi: 10.2967/jnumed.109.069732
21. Nakamuta JS, Danoviz ME, Marques FLN, Dos Santos L, Becker C, Goncalves GA, et al. Cell therapy attenuates cardiac dysfunction post myocardial infarction: effect of timing, routes of injection and a fibrin scaffold. PLoS ONE. (2009) 4:e6005. doi: 10.1371/journal.pone.0006005
22. Swijnenburg R-J, Govaert JA, Van Der Bogt KEA, Pearl JI, Huang M, Stein W, et al. Timing of bone marrow cell delivery has minimal effects on cell viability and cardiac recovery after myocardial infarction. Circulation. (2010) 3:77–85. doi: 10.1161/CIRCIMAGING.109.872085
24. Kadner K, Dobner S, Franz T, Bezuidenhout D, Sirry MS, Zilla P, et al. The beneficial effects of deferred delivery on the efficiency of hydrogel therapy post myocardial infarction. Biomaterials. (2012) 33:2060–6. doi: 10.1016/j.biomaterials.2011.11.031
25. Bastings MMC, Koudstaal S, Kieltyka RE, Nakano Y, Pape ACH, Feyen DAM, et al. A fast pH-switchable and self-healing supramolecular hydrogel carrier for guided, local catheter injection in the infarcted myocardium. Adv Healthc Mater. (2014) 3:70–8. doi: 10.1002/adhm.201300076
26. Wise P, Davies NH, Sirry MS, Kortsmit J, Dubuis L, Chai C-K, et al. Excessive volume of hydrogel injectates may compromise the efficacy for the treatment of acute myocardial infarction. Int J Numer Method Biomed Eng. (2016) 32:e02772. doi: 10.1002/cnm.2772
27. Dobner S, Bezuidenhout D, Govender P, Zilla P, Davies N. A synthetic non-degradable polyethylene glycol hydrogel retards adverse post-infarct left ventricular remodeling. J Card Fail. (2009) 15:629–36. doi: 10.1016/j.cardfail.2009.03.003
28. Ciuffreda MC, Malpasso G, Chokoza C, Goetsch KP, Mura M, Pisano F, et al. Synthetic extracellular matrix mimic hydrogel improves efficacy of mesenchymal stromal cell therapy for ischemic cardiomyopathy. Acta Biomater. (2018) 70:71–83. doi: 10.1016/j.actbio.2018.01.005
29. Wall ST, Walker JC, Healy KE, Ratcliffe MB, Guccione JM. Theoretical impact of the injection of material into the myocardium. A finite element model simulation. Circulation. (2006) 114:2627–35. doi: 10.1161/CIRCULATIONAHA.106.657270
30. Wenk JF, Eslami P, Zhang Z, Xu C, Kuhl E, Gorman JH, et al. A novel method for quantifying the in-vivo mechanical effect of material injected into a myocardial infarction. Ann Thorac Surg. (2011) 92:935–41. doi: 10.1016/j.athoracsur.2011.04.089
31. Kortsmit J, Davies NH, Miller R, Macadangdang JR, Zilla P, Franz T. The effect of hydrogel injection on cardiac function and myocardial mechanics in a computational post-infarction model. Comput Methods Biomech Biomed Engin. (2013) 16:1185–95. doi: 10.1080/10255842.2012.656611
32. Kortsmit J, Davies NH, Miller R, Zilla P, Franz T. Computational predictions of improved of wall mechanics and function of the infarcted left ventricle at early and late remodelling stages : comparison of layered and bulk hydrogel injectates. Adv Biomech Appl. (2013) 1:41–55. doi: 10.12989/aba.2013.1.1.041
33. Miller R, Davies NH, Kortsmit J, Zilla P. Outcomes of myocardial infarction hydrogel injection therapy in the human left ventricle dependent on injectate distribution. Int J Numer Methods Biomed Eng. (2013) 29:870–84. doi: 10.1002/cnm.2551
34. Sirry MS, Davies NH, Kadner K, Dubuis L, Saleh MG, Meintjes EM. Micro-structurally detailed model of a therapeutic hydrogel injectate in a rat biventricular cardiac geometry for computational simulations. Comput Methods Biomech Biomed Engin. (2013) 18:325–31. doi: 10.1080/10255842.2013.793765
35. Kichula ET, Wang H, Dorsey SM, Szczesny SE, Elliot DM, Burdick JA, et al. Experimental and computational investigation of altered mechanical properties in myocardium after hydrogel injection. Ann Biomed Eng. (2014) 42:1546–56. doi: 10.1007/s10439-013-0937-9
36. Lee LC, Wall ST, Genet M, Hinson A, Guccione JM. Bioinjection treatment: effects of post-injection residual stress on left ventricular wall stress. J Biomech. (2014) 47:3115–9. doi: 10.1016/j.jbiomech.2014.06.026
37. Sack KL, Davies NH, Guccione JM, Franz T, Franz T. Personalised computational cardiology : patient-specific modelling in cardiac mechanics and biomaterial injection therapies for myocardial infarction. Heart Fail Rev. (2016) 21:815–26. doi: 10.1007/s10741-016-9528-9
38. Daviran M, Longwill SM, Casella JF, Schultz KM. Rheological characterization of dynamic remodeling of the pericellular region by human mesenchymal stem cell-secreted enzymes in well-defined synthetic hydrogel scaffolds. Soft Matter. (2018) 14:3078–89. doi: 10.1039/C8SM00408K
39. Poudel AJ, He F, Huang L, Xiao L, Yang G. Supramolecular hydrogels based on poly (ethylene glycol)-poly (lactic acid) block copolymer micelles and α -cyclodextrin for potential injectable drug delivery system. Carbohydr Polym. (2018) 194:69–79. doi: 10.1016/j.carbpol.2018.04.035
40. Naghizadeh Z, Karkhaneh A, Khojasteh A. Self-crosslinking effect of chitosan and gelatin on alginate based hydrogels: injectable in situ forming scaffolds. Mater Sci Eng C. (2018) 89:256–64. doi: 10.1016/j.msec.2018.04.018
41. Holmes R, Yang XB, Dunne A, Florea L, Wood D, Tronci G. Thiol-Ene Photo-Click Collagen-PEG hydrogels: impact of water-soluble photoinitiators on cell viability, gelation kinetics and rheological properties. Polymers. (2017) 9:226. doi: 10.3390/polym9060226
43. Shulyak I V, Grushova EI, Semenchenko AM. Rheological properties of aqueous solutions of polyethylene glycols with various molecular weights. Russ J Phys Chem A. (2011) 85:485–8. doi: 10.1134/S0036024411030265
44. Baldursdottir S, Tian F, Santacruz BG, Jørgensen AC. The influence of thermal history on the physical behavior of poly (ethylene glycol) (PEG). Pharm Dev Technol. (2012) 17:195–203. doi: 10.3109/10837450.2010.531733
45. Shulyak IV, Grushova EI. Rheological properties of water – salt solutions of polyethylene glycols of various molecular masses in the range of 293.15 to 323.15 K. Russ J Phys Chem A. (2013) 87:453–6. doi: 10.1134/S0036024413030291
46. Azri A, Privat M, Grohens Y, Aubry T. Linear rheological properties of low molecular weight polyethylene glycol solutions. J Colloid Interface Sci. (2013) 393:104–8. doi: 10.1016/j.jcis.2012.10.059
50. Alonso C, Pries AR, Kiesslich O, Lerche D, Gaehtgens P. Transient rheological behavior of blood in low-shear tube flow: velocity profiles and effective viscosity. Am J Physiol. (1995) 268:H25–32. doi: 10.1152/ajpheart.1995.268.1.H25
55. Sherwood JM, Kaliviotis E, Dusting J, Balabani S. Hematocrit, viscosity and velocity distributions of aggregating and non-aggregating blood in a bifurcating microchannel. Biomech Model Mechanobiol. (2014) 13:259–73. doi: 10.1007/s10237-012-0449-9
58. Sands GB, Gerneke DA, Hooks DA, Green CR, Smaill BH, Legrice IJ. Automated imaging of extended tissue volumes using confocal microscopy. Microsc Res Tech. (2005) 67:227–39. doi: 10.1002/jemt.20200
60. Merck. Polyethylene glycol 20000 SDS. Merck Saf Data Sheet. (2018). Available online at: http://www.merckmillipore.com/ZA/en/product/msds/MDA_CHEM-818897?Origin=PDP (accessed April 30, 2019).
61. Van den Akker F, Feyen DAM, Van Den Hoogen P, Van Laake LW, Van Eeuwijk ECM, Hoefer I, et al. Intramyocardial stem cell injection: Go(ne) with the flow. Eur Heart J. (2017) 38:184–6. doi: 10.1093/eurheartj/ehw056
Keywords: myocardial infarction, injectate therapy, polyethylene glycol hydrogel, particle image velocimetry, computational fluid dynamics, retention
Citation: Ngoepe M, Passos A, Balabani S, King J, Lynn A, Moodley J, Swanson L, Bezuidenhout D, Davies NH and Franz T (2019) A Preliminary Computational Investigation Into the Flow of PEG in Rat Myocardial Tissue for Regenerative Therapy. Front. Cardiovasc. Med. 6:104. doi: 10.3389/fcvm.2019.00104
Received: 12 May 2019; Accepted: 16 July 2019;
Published: 07 August 2019.
Edited by:Rajesh Katare, University of Otago, New Zealand
Reviewed by:Changyong Wang, Beijing Institute of Basic Medical Sciences, China
Federico Quaini, University of Parma, Italy
Copyright © 2019 Ngoepe, Passos, Balabani, King, Lynn, Moodley, Swanson, Bezuidenhout, Davies and Franz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Malebogo Ngoepe, email@example.com
†These authors have contributed equally to this work