Mathematical modeling and experimental testing of three bioreactor configurations based on windkessel models.

This paper presents an experimental study of three bioreactor configurations. The bioreactor is intended to be used for the development of tissue-engineered heart valve substitutes. Therefore it must be able to reproduce physiological flow and pressure waveforms accurately. A detailed analysis of three bioreactor arrangements is presented using mathematical models based on the windkessel (WK) approach. First, a review of the many applications of this approach in medical studies enhances its fundamental nature and its usefulness. Then the models are developed with reference to the actual components of the bioreactor. This study emphasizes different conflicting issues arising in the design process of a bioreactor for biomedical purposes, where an optimization process is essential to reach a compromise satisfying all conditions. Two important aspects are the need for a simple system providing ease of use and long-term sterility, opposed to the need for an advanced (thus more complex) architecture capable of a more accurate reproduction of the physiological environment. Three classic WK architectures are analyzed, and experimental results enhance the advantages and limitations of each one.

Introduction

Bioreactors are widely used in research involving vascular and heart valve substitutes, whether these substitutes are made of biomaterials or constructed using tissue engineering techniques. They are necessary to reproduce physiological conditions, either for the observation of the performance of these substitutes or for the elaboration of new techniques yielding these constructs. A bioreactor designed for this purpose should allow the reproduction of blood flow and arterial pressure, as well as other relevant biochemical and physical parameters such as temperature, pH, CO2, and O2 concentrations. Long-term sterility is also an important issue. Difficulties associated with the complexity of the parameters to simulate render an exact reproduction of this environment difficult. The foremost difficulty arises from the dynamic pulsating nature of arterial blood flow and pressure. The strongly nonlinear relation existing between these two quantities is a direct consequence of the nonlinearity of the involved physical phenomena, notably the interaction between the fluid and the visco-elastic tissues. The global physics of the fluid flow produced by a bioreactor is often studied using mathematical models consisting of different assemblies of basic components, and are referred to as windkessel (WK) or lumped models. The basic components used to represent the actual bioreactor and translate it into a mathematical model are a flow source, flow resistances, elastic accumulators, and inertial components representing the mass inertia of blood in arteries. Each one of these elements represents a physical component of the bioreactor, which in turn simulates an actual property of the cardiovascular environment.

Bioreactors generally feature a flow source producing pulsating flow rates. A resistance (this parameter is noted R in WK models) mimics the blood flow restrictions caused by arteries, capillaries, and veins. It is responsible for systemic pressure increase when flow is forced into the circulatory system by the pump. The elastic accumulator is obtained by the insertion, in the bioreactor circuit, of a closed volume containing a certain amount of air. During systole, the air is compressed in the accumulator because the fluid is prevented from circulating freely by the flow resistance. The flow accumulator thus simulates the physiological compliance property (C parameter in WK models) associated with tissue proclivity to expand under a blood pressure increase. Many studies have reported the importance of an accurate reproduction of compliant effects in bioreactors,[1,2] especially in tissue engineering applications[3-6] where strong evidence of a significant relation between the dynamic frequency spectrum of the pulsating flow and tissue compliance has been observed.

The work presented in this paper has been carried out in the design process of a new bioreactor dedicated specifically to tissue engineering of heart valves. One of the objectives underlying the project was an optimal simultaneous reproduction of the dynamic waveforms of arterial pressure and blood flow rate. The aim of this article is not to present the bioreactor (which we have reported[7]), but to present experimental data allowing a physical comparison between three much-covered WK architectures. No such systematic, experiment-based comparison is available in the literature, although a very large number of articles report WK-related studies.

Windkessel models in the literature

Previous bioreactors producing physiological flow and pressure waveforms were either explicitly based on WK models, such as those designed by Dumont et al.,[8] Hildebrand et al.,[9] and Knierbein et al.,[10] or entail the appropriate elements structured in an adequate manner to simulate the hemodynamic environment, like those developed by Warnock et al.,[11] Konduri et al.,[12,13] and Narita et al.[14] Generally, macro-physiological properties of the cardiovascular system, such as the aortic and arterial compliance, resistance, and inertial effects, are modeled using few components (from two to four elements). It can be noticed in several articles presenting mathematical analyses of different arrangements that the complexity of the differential equation describing the system rapidly increases with the number of elements constituting the model. Some papers propose a dynamic response/transfer function type analysis,[15,16] while others directly solve the differential equation by numerical methods.[17,18] Stergiopulos et al.[19-21] studied the differences between three- and four-component architectures and compared their theoretical responses to physiological data, but not to experimental data obtained with actual bioreactors. More complex models of the arterial tree have also been proposed and combined to numerical resolution algorithms, providing another comprehensive tool for the study of cardiovascular behavior. Conlon et al.[22] published a model composed of 14 elements whereas another model, developed by Liang et al.,[23] included more than 150 elements. The usefulness of the WK approach is also enhanced in numerous experimental studies addressing physiological data analysis of human subjects. In particular, it has been demonstrated that frequency content analysis of arterial pressure waveforms and WK modeling could advantageously support pathology identification techniques, notably for hypertension.[7,24-28]

The literature review of papers presenting the WK approach permits one to appreciate the depth of the concept as well as its ability to model the dynamics of the cardiovascular system, especially in cases of models that include more than two components. Although many papers outline the assets of WK architectures more complex than the classic two-component model, the vast majority of bioreactors found in publications are actually developed according to that basic configuration. Some papers present a detailed theoretical study, sometimes with comparisons to physiological data.[19,20,22] Others present sophisticated bioreactor configurations[8,9,11,14] but without a systematic theoretical analysis to support and justify the design. In the present article, we propose a bridge between WK modeling, theoretical study, and the actual construction of corresponding bioreactors. We will present a complete and systematic analysis of the three most common WK architectures: RC, RRC, and RRCL. In these designations, R and C represent the resistive and compliant components, while L represents inertia. A recent article by Westerhof et al.[29] presents a very interesting review of theoretical aspects related to these three WK models. In the work presented here, the differential equation will be established and solved numerically for each WK architecture. Theoretical results will be compared to experimental data obtained from three actual bioreactors that were built according to these architectures. This will establish the benefit of bioreactor architecture complexity in providing an accurate simultaneous reproduction of both the flow and the pressure waveforms.

Materials and Methods

RC bioreactor

Schematic representations of the first bioreactor configuration (RC model) are presented in Figure 1. The pump drives the flow in the circuit, composed of a compliance chamber C followed by a resistance R and a return to the pump. C and R represent the properties of the arterial tree considered in its entirety. The test section is located immediately after the pump, so that the conditions felt by the heart valve are a flow rate q and a pressure P. Indices indicate in which component a physical parameter is evaluated. A corresponding electrical circuit is also presented in Figure 1, as this is the simplified representation generally encountered in the literature related to WK models, and also because of the physical analogy between the two systems. Tissue compliance, simulated by the compliance chamber (and analog to an electric capacitance) is defined as the ratio between the difference in volume induced by a pressure variation and that pressure variation:[10,16] where V is the volume and PC the pressure in the compliance chamber. Deriving numerator and denominator with respect to time, a direct relation is obtained between the flow and the pressure gradient: or

Figure 1

(A) RC bioreactor configuration; (B) corresponding electric circuit.

Assuming a negligible pressure drop between the pump and the compliance chamber (very short length of tube linking these two elements), we have

The analysis is based on the assumption that the flow is laminar (which is generally true for physiological flows), and the relation between q and the pressure difference ΔP across the resistance can be expressed as: where cst is a constant. Defining the hydraulic resistance as the ratio of the pressure increase to the imposed flow, we have

Assuming that the pressure difference between the exit of the resistance and the entrance of the pump is negligible, we have: Then

Noting that q = q + q, the following differential equation is obtained for the RC bioreactor:

An actual bioreactor based on this two-component RC configuration was built.[7] The flow is produced by a custom design diaphragm pump, driven by a computer-controlled proportional servo-valve (Festo MPYE-5-⅛, Mississauga, ON, USA). The pump allows a pulsating frequency range of 60–180 bpm and a stroke volume up to 160 mL. From the exit of the pump, the fluid flows through a flexible 25-mm internal diameter tube (Tygon, US Plastic Corp., Lima, OH, USA), onto which is mounted a noninvasive ultrasonic flow sensor probe (Transonics Systems Inc., Ithaca, NY, H20XL). Immediately downstream, a valve holder maintains the aortic valve while exposing it to a surrounding geometry reproducing the complex tridimensional sinuses forming the aortic root, thus providing an accurate reproduction of physiological flow dynamics. Fluid pressure upstream and downstream of the valve is monitored with 0–7 psi dynamic pressure sensors (Argon Medical Devices, Athens, TX, USA). A compliance chamber (custom-made glass reservoir) is located after the valve holder, and followed by a variable resistance. The compliance value can be adjusted between 0 and 1 mL/mmHg, while the resistance range is 0.5–1.2 mmHg.s.mL−1. All components and sensors mentioned in this section are interfaced to a USB National Instruments I/O board, with control and data acquisition tasks performed through a custom-built LabView application program.

The difference between this arrangement (Figure 2) and the previous one is that a second resistance RC has been added between the pump and the compliance chamber. The resistance Rp represents the peripheral resistance of the complete arterial tree. The physical meaning of the added resistance RC is more subtle. Although RC presents the units of a resistance, it corresponds to the characteristic (or total) impedance[16] of the aorta, a combination of its local compliance and the inertance of the blood volume. According to the continuity equation, the relation between the different flow rates is:

Figure 2

(A) RRC bioreactor configuration; (B) corresponding electric circuit.

Flow qRc can be expressed by the ratio of the pressure difference across the resistance RC, and the resistance RC itself:

Replacing terms qC, qRc, and qRp in equation [10] by their corresponding expressions yields the following differential equation:

Assuming that the pressure difference between the exit of the resistance and the entrance of the pump is negligible, that is: the following differential equation is obtained for the RRC bioreactor:

To obtain the RRC arrangement with the actual bioreactor, a second resistance was added between the valve holder and the compliance chamber, while other components remained unchanged.

Schematic representations of the RRCL bioreactor are presented in Figure 3. Here a tube with a certain length l and a cross-section A is added in parallel to the resistance RC. The flow restriction in this tube is low compared to that of the resistance RC. However, the fluid mass it contains presents inertia opposed to large accelerations. As outlined by Stergiopulos et al.,[20] the frequency response of this arrangement is more likely to produce a physiological pressure waveform corresponding to the one normally observed in vivo; that is, the natural response to a physiological flow rate input. For low frequencies (near steady flow conditions), the fluid flows freely in the tube as its resistance is lower than that of RC. On the other hand, higher frequency components of the flow rate cannot circulate freely because fluid inertia opposes sudden accelerations: the higher frequencies therefore flow through the resistance RC and are attenuated in amplitude. The added tubing thus acts like an electric inductance L. According to the continuity equation, the relation between the different flow rates is:

Figure 3

(A) RRCL bioreactor configuration; (B) corresponding electric circuit.

Flow qL is defined using the relation existing between pressure and flow modulated by inertial effects (inertance L) of the moving mass of fluid:

This relation between flow and pressure is a direct result from the fact that a pressure difference ΔPL is required to produce acceleration. Considering inertial effects as widely predominant with respect to viscous friction at the wall, we can express this relation with Newton’s second law: which can be easily expressed in terms of the flow parameters: where ρ is the fluid’s density, and the tube’s cross section A is assumed to be constant over time. Rearranging equation [18], we obtain an expression for the flow rate qL:

Replacing the terms q and qL in equation [15] by their corresponding expressions [equations 3, 11, and 16] yields the following differential equation: and noting that the following differential equation is obtained for the RRCL bioreactor:

To obtain the RRCL arrangement with the actual bioreactor, flow dividers were inserted before and after the resistance RC, and a certain length of tubing was connected in parallel to that resistance.

To test the three configurations presented in the previous section, reference waveforms corresponding to clinical measurements of arterial flow rate and pressure performed simultaneously on a healthy individual[12] were used. Calculations with these models were widely used in the design process, with a constant retroaction between theory and design to dimension bioreactor components properly. The first step was to perform, for each model, a parametric study enhancing the influence of the R, C, and L parameters on the system response. A Matlab program was developed to solve the differential equations [9], [14] and [22] corresponding to the three models and perform these parametric studies. For this the flow rate (dotted line in Figure 4A) was used as the input of the system, and the corresponding pressure response was calculated as a function of R, C, and L. Optimal parameter values were obtained by implementing the conjugate gradients method in the differential equation solver program to determine, from a theoretical point of view, which parameters were more likely to produce the reference pressure waveform when the reference flow rate waveform was used as the input. Then the dimensions of the bioreactor were calculated from these results to physically achieve these optimal parameter values. The bioreactor was constructed and the three arrangements were tested and compared, first to the predictions of their corresponding mathematical model, then to each other.

Figure 4

Experimental results for the RC configuration.

To compare the waveforms and obtain an evaluation of the correlation between two curves, a coefficient of correlation r was calculated: where P represents discrete experimental data, P is the mean value of the experimental data, P represents discrete reference data, and P is the mean value of the reference data. The coefficient of correlation r ranges from 0–1, the upper limit corresponding to a perfect fit.

Results and Discussion

Figure 4 presents experimental measurements of the flow rate and pressure waveforms produced by the RC bioreactor. Water was used for these experiments. It can be observed that the shape of the measured flow rate (Figure 4A) is very similar to the physiological reference waveform, and the coefficient of correlation between the two curves is 0.976. The pressure waveform measured in the bioreactor (Figure 4B) is different from the reference pressure, but nevertheless presenting the same general shape. The differences between the waveforms are a certain phase lag as well as a higher amplitude. The delay is related to the fact that air must be compressed in the compliant element before the pressure increases, while the higher amplitude can be attributed to a slightly suboptimal compliance value. The coefficient of correlation between the measured pressure and the reference is 0.335, owed mainly to the delay rather than to the shape of the waveform. The theoretical RC pressure response, calculated by optimizing the model to obtain experimental R and C parameters (R=1.16 mmHg.s.mL−1 and C=0.91 ml.mmHg−1), is in good agreement with the experimental pressure waveform (coefficient of correlation = 0.960). The higher frequency oscillations of the pressure waveform at minimal and maximal values may be attributed to vibrations of the flexible components of the system (tubing) at the opening and closing of the valve leaflets, and by free surface movement in the compliance chamber. These results show that appropriately calibrated resistance and compliance components produce a good first approximation of the pressure response, and also that the WK model describes the global physics of the actual bioreactor well (excellent correlation between model and measurements) and constitutes an appropriate approach to support its design and dimensioning.

In the RRC configuration, the addition of a second resistance to the circuit before the compliance chamber eliminates the phase shift: a direct and simultaneous pressure increase is produced (proportionally) by the passage of the flow through the resistance. This can be appreciated in Figure 5, where the experimental pressure curve can be seen in phase with the physiological reference waveform. The coefficient of correlation is improved to a value of 0.766. Excellent agreement can be observed between the theoretical RRC model and the experimental results (coefficient of correlation = 0.947).

Figure 5

Experimental results for the RRC configuration.

The addition of a fourth component further improves the natural response of the bioreactor (Figure 6). The value of the coefficient of correlation is increased to 0.816, and the overall pressure waveform presents a more accurate reproduction of the shape characteristics of the physiological reference. In particular, the reproduction of the dicrotic pulse at the end of systole can be observed. The inertance L added in parallel to the resistance RC located upstream from the compliance chamber complicates the dynamic interaction between RC, L, and C. It was found that a significant length of tubing was required to achieve the optimal response of the system (L ≈ 2.5 m). At the beginning of systole, the fluid must be set in motion. If it flows through the resistance RC, the pressure increases proportionally to the flow. If it flows through the inertance L, mass inertia also produces a pressure increase. The balance between the amount of fluid flowing through RC and L depends on the relative magnitude of the impedance of both branches (resistance and inertance). A large resistance value will force the fluid to flow through the inertance branch if its length L is such that this path is easier. If L is very long and the mass is difficult to set in motion, the fluid may prefer to flow through the resistance. This balance between the two paths changes dynamically because when the fluid starts to flow into the inertance branch, its resistance to motion decreases and requires a lower pressure to maintain the flow.

Figure 6

Experimental results for the RRCL configuration.

For all combinations of RC and L, it was observed experimentally that the initial phase of systole produced a pressure increase in good agreement with the physiological reference. Then a combination of the two following phenomena occurred: when the fluid gradually gained velocity in the inductance branch, the flow resistance decreased and the pressure decreased accordingly as this path offered less restriction; then a flow inversion occurred in the RC branch and fed the inductance branch. The momentum was simply transferred quite freely from one branch to another and the consequence on the pressure was a truncation of the systole pressure peak. It was observed that low resistance values enhance backflow occurrence and truncation of the systole pressure peak. To avoid this truncation during the systolic phase, it was found that a high RC resistance value was required. In addition, to avoid fluid bypass of the RC branch, a large inductance value was required. Consequently, high settings of both RC and L ensured that sufficient flow resistance was present in both branches, resulting in an adequate pressure increase over the whole systolic phase. These high settings of RC and L also meant that the flow inverting phenomena in the RC branch was less pronounced in amplitude as well as more delayed in time, but nevertheless present. In fact it occurred later, in the decreasing part of the systolic waveforms, and yielded an additional momentum pulse to an otherwise decreasing waveform, thus reproducing the dicrotic pulse.

For the third configuration, the correlation between the theoretical model and the experimental results was good but lower than for the two previous architectures (coefficient of correlation = 0.741). This can be attributed to the fact that the differences between the ideal and actual components are more important; namely, because of the connections geometry and also because of the multiple-nature of RC and L. In reality L is not purely inductive as it also offers a certain resistance, and correspondingly RC is not purely resistive. This emphasizes certain limitations of WK models in providing accurate information about actual physical systems.

Even considering the fact that the four-component architecture improves the natural, passive response of the bioreactor, we finally decided that future work would be based on the simpler three-component configuration for two main reasons. First, our aim is to reach the maximum level of simplicity, because fewer components mean easier assembly and more efficient sterilization operations, which in turn should help long-term sterility. This is a major issue for a bioreactor designed specifically for the development of tissue engineering construction techniques, which may require up to several weeks of continuous operation. Second, the modification of flow conditions with a four-component bioreactor requires changing hardware components, which is obviously more complicated (if not completely impossible without interrupting flow circulation) than just programing a simpler setup that could eventually be equipped with controllable components.

Conclusions

In this paper we presented a detailed study of the theoretical response of three WK configurations. We then compared these results to experimental data obtained with three actual bioreactors constructed according to these architectures. These results showed that: i) WK models are good tools to support the design of experimental bioreactors, and ii) the addition of a third and a fourth component to the bioreactor greatly improved the passive reproduction of the pressure waveform to a physiological flow rate.