Ultrasound in Med. & Biol., Vol. 41, No. 6, pp. 1766–1778, 2015 Copyright
2015 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/$ - see front matter

http://dx.doi.org/10.1016/j.ultrasmedbio.2014.12.015

d

Original Contribution THEORETICAL EVALUATION OF THE ACOUSTIC FIELD IN AN ULTRASONIC BIOREACTOR TOBIAS M. LOUW, ANURADHA SUBRAMANIAN, and HENDRIK J. VILJOEN Department of Chemical and Biomolecular Engineering, University of Nebraska—Lincoln, Lincoln, Nebraska, USA (Received 18 July 2014; revised 11 December 2014; in final form 15 December 2014)

Abstract—Ultrasound-assisted bioreactors that provide mechanical conditioning to cells have broad applicability in tissue engineering, but biological experiments with ultrasound are very sensitive to environmental conditions. A mathematical model was developed to complement experimental measurements, as well as to describe ultrasonic fields existing in regions where measurements are impossible, specifically, within microporous tissue engineering scaffolds. The model uniquely combines Biot theory to predict the ultrasonic field in the scaffold with an electromechanical transducer model to couple the mechanical stimulation experienced by cells to the external electrical input. In the specific example examined here, cells immobilized on scaffolds are subjected to different forms of ultrasonic stimulation due to the formation of standing wave fields and vertical high-pressure bands. The model confirms the sensitivity of the supplied acoustic power to the liquid level in sonobioreactors and identifies the input electrical impedance as a method of detecting resonance effects. (E-mail: [email protected])
2015 World Federation for Ultrasound in Medicine & Biology. Key Words: Ultrasound, Bio-acoustics, Mechanotransduction, Tissue engineering, Finite volume method, Spectral method, Biot theory.

affects cells is not well understood (Dalecki 2004; Fowlkes 2008; Nelson et al. 2009). Sonobioreactors and the ability to conduct reproducible in vitro experiments are central to tissue engineering and medical ultrasound research in general. The effects of ultrasound are usually investigated by stimulating a group of cells and measuring the average biological response. It is important to ensure that all cells in the treatment group receive similar acoustic simulation to accurately correlate cellular responses to ultrasonic stimulation. Furthermore, acoustic regimens must be controlled very precisely to avoid negative bio-effects induced by cavitation or thermal events (Karande et al. 2005; Kinoshita and Hynynen 2007; Mitragotri 2005). Ultrasound wave propagation in routinely used in vitro setups have been analyzed by using either a combination of simulations and hydrophone measurements (Hensel et al. 2011) or a combination of pulse-echo ultrasound, laser Doppler vibrometry and Schlieren imaging (Leskinen and Hynynen 2012). These studies highlighted the sensitivity of experiments in ultrasonic bioreactors to external disturbances, with frequency effects creating uncertainties of up to 700%, and concluded that detailed analyses are required to control experimental variations (Leskinen and Hynynen 2012).

INTRODUCTION Ultrasound is frequently used in therapeutic applications such as wound healing, targeted drug delivery and lithotripsy (Dalecki 2004; Naito et al. 2010; Rubin et al. 2001), as well as in in vitro cultures (Hsu et al. 2006a, 2006b; Nishikori et al. 2002; Parvizi et al. 1999; Zhang et al. 2003). Ultrasonic bioreactors (also known as sonobioreactors) (Fig. 1) can provide the necessary mechanical stimulation required to reproduce physiologic conditions and promote synthesis of biological tissue substitutes (Hsu et al. 2007; Whitney and Lamb 2012; Zhang et al. 2003). Cells are most often stimulated by pulsed ultrasound, but recent work has indicated that the stimulation of in vitro chondrocyte cultures by continuous ultrasound has a positive influence on several factors, including cell proliferation, viability and gene expression of select chondrocytic and load inducible-markers (Hasanova et al. 2011; Louw et al. 2013; Noriega et al. 2013; Whitney and Lamb 2012). Despite the many applications of ultrasound in biology, the exact mechanism whereby ultrasound Address correspondence to: Hendrik J. Viljoen, Department of Chemical and Biomolecular Engineering, University of Nebraska— Lincoln, 207 Othmer Hall, 820 North 16th Street, Lincoln, NE 68588, USA. E-mail: [email protected] 1766

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

1767

Our goals in this work were to (i) develop a framework for applying an appropriate mathematical model with the ability to predict the propagation of ultrasound through liquids, elastic solids and porous media, to a sonobioreactor; (ii) identify and validate appropriate numerical methods for solving the mathematical model; (iii) couple the model-predicted ultrasonic field to the measureable transducer electrical input; (iv) use the mathematical model to identify possible disturbances that may cause experimental error by affecting the ultrasonic field; and (v) identify potential measurement and control variables as well as design criteria that can be used to ensure repeatability of sonobioreactor experiments. An accurate and efficient mathematical model can be used to design experiments aimed at discovering the biological effects of ultrasound and optimize bioreactors for the eventual production of synthetic tissues and other high-value biological products. Fig. 1. Typical sonobioreactor setup. (a) The entire setup, excluding electrical components, is kept in an incubator with temperature and oxygen controlled. (b) Cell seeded scaffolds are placed in the wells of a tissue culture plate; each well contains growth medium. Tissue culture plates can contain 6, 12 or 24 wells. (c) The tissue culture plate is partially submerged in an aquarium so that no air bubbles exist below the tissue culture wells, while the culture medium is kept isolated from the aquarium. The aquarium acoustically couples the wells to the ultrasonic transducers. (d) An ultrasonic immersion transducer is placed below each of the tissue culture wells to provide mechanical stimulation (bottom-up sonication). Alternatively, the transducers may be placed directly above the tissue culture well, partially immersed in the culture medium, for top-down sonication. (e) A function generator (FG) and current splitter provide the voltage input to each ultrasound transducer. The input is sinusoidal with variable frequency and voltage amplitude.

A mathematical model is one such analysis method. Such a model would be especially useful in predicting ultrasonic effects within regions inaccessible to measurement, such as the interior of tissue engineering scaffolds. The original theory for the propagation of sound through porous media was developed by M.A. Biot (1962), and is valid whenever the acoustic wavelength is much larger than the characteristic dimensions of the porous structure under investigation, as is the case here. Furthermore, the bioreactor is driven by a real transducer with a complex beam shape that must be accounted for. The power delivered to the bioreactor by the transducer is dependent on the electrical input as well as the transducer impedance. It is necessary to couple the bioreactor model to an electromechanical transducer model to quantitatively predict the ultrasonic field in the bioreactor and scaffold. A unique combination of mathematical models, including the effects of porous media, real transducers and electrical inputs, is necessary to model an ultrasonic bioreactor.

METHODS Mathematical modeling An axisymmetric mathematical model suffices for the bioreactor setup illustrated in Figure 2(a). The bioreactor consists of fluids (water in aquarium, culture media), solids (tissue culture well) and poro-elastic media (scaffolds). The ultrasonic transducer acts as the acoustic source and fixes the fluid velocity v 5 ½vr ; vz  in the z-direction such that vz ðr; z; tÞ 5 v0 ðtÞ at the transducer face. The velocity amplitude v0 depends on the electromechanical coupling occurring between the supplied voltage source and the transducer face by means of a piezoelectric slab, as predicted by the Krimholtz–Leedom–Matthaei (KLM) transducer model (Krimholtz et al. 1970). Two numerical methods are investigated and compared in terms of accuracy and computational efficiency. The transfer matrix/angular spectrum approach (TM/ASA) is a semi-analytical method acting under the fundamental assumption that acoustic wave reflections in the radial direction can be neglected (Louw 2013). Mathematically, the TM/ASA represents the ultrasonic transducer as being set in an infinite, rigid baffle and approximates every bioreactor component (scaffold, tissue culture well, etc.) as infinitely wide, hence no radial reflections (see Fig. 2b). In practice, the assumption is justified when each acoustic component has a diameter greater than that of the ultrasonic beam; this postulate is evaluated by comparing TM/ ASA results to the output of a second numerical method: the finite volume method (FVM) (Saenger et al. 2000; Virieux 1986). The FVM is a full numerical method that accounts for the complex geometry of the entire system, but has much higher computational cost than the TM/ASA. By comparing

1768

Ultrasound in Medicine and Biology

Volume 41, Number 6, 2015

Fig. 2. Diagram of the sonobioreactor, with boundary conditions. (a) Scaffolds are placed in a cylindrical polystyrene tissue culture well. Axial symmetry is assumed in the mathematical model. Different setups include with/without scaffold, with/without hydrophone and transducer placement below/above tissue culture well. Boundary conditions are indicated. Zh and Znr are the characteristic acoustic impedances of the hydrophone and the medium at the nonreflecting boundary, respectively. (b) The transfer matrix/angular spectrum approach (TM/ASA) approximates each acoustic layer as being infinitely wide; hence all radial reflections and boundary conditions are neglected.

the TM/ASA results with those generated by the FVM, the assumption of negligible radial reflection underlying the TM/ASA can be validated. The TM/ ASA and FVM are merely two different methods used to solve the equations describing the ultrasonic field in the sonobioreactor; the next section describes these governing equations.




"

vt sðsÞ vt sðfÞ

Ultrasound in inviscid fluids, linear elastic solids and poro-elastic media The momentum balance and infinitesimal stressstrain relationships that govern wave propagation in linear elastic solids are calculated as r

vv 5 V$s vt

vs vU 5 ll ðV$vÞI12G vt vt

(1)

(2)

where v is the velocity vector, and s is the stress tensor (Allard and Atalla 2009). In eqns (1) and (2), r is the density, ll is the first Lame parameter, G is the shear modulus, U is the small strain tensor and I is the 333 identity matrix. Substituting s 5 --pI (pressure), ll 5 K (bulk modulus) and G 5 0 (no shear stress) in eqns (1) and (2) yields the linearized momentum balance and continuity equation for inviscid fluids. The corresponding equations for poro-elastic solids are (Biot 1962)

  ðsÞ    r12 vt v V$sðsÞ 5I 5 r22 vt vðfÞ V$sðfÞ

r11 r12

#

 5

ðP22NÞ

Q

Q

R

5I12N

"

vt U ðsÞ

" #

V$vðsÞ V$vðfÞ

(3)

#! (4)

0

Superscripts ðsÞ and ðfÞ refer to the solid and fluid phases, respectively, 5 is the Kronecker product and 0 is a 333 matrix of zeroes. Assuming harmonic time dependence (because of the application of continuous ultrasound), the time derivative operator vt can be replaced by iu. The densities rij and parameters P; Q; R and N are material properties of the poro-elastic medium, as explained in the Supplementary Material (see online version at http://dx.doi.org/10.1016/j.ultrasmedbio.2014.12.015) (Allard and Atalla 2009). Coupling conditions at the interface between different media ensure continuity of stress, fluid flow and solid displacement. The exact equations describing the coupling conditions are given in the Supplementary Material. With the application of boundary conditions (see Fig. 2a), the governing equations (eqns 1–4) and coupling conditions form a fully determinate system of partial differential equations. Two different methods are considered to solve the set of equations.

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

Transfer matrix/angular spectrum approach The TM/ASA reduces the system of partial differential equations to a set of linear algebraic equations by representing the ultrasonic standing wave field as the sum of infinitely many plane waves. A real transducer in a fluid medium generates an intricate pressure field pðr; zÞ, described by the Rayleigh–Sommerfeld integral (Kinsler et al. 1999). The first step in the TM/ASA is to apply a Fourier transform ℱ f$g to the pressure field predicted by the Rayleigh–Sommerfeld integral at the transducer face, z 5 0. The acoustic field at z 5 0 can be viewed as the superposition of infinitely many plane waves, identified by the continuous variable x˛ð2N; NÞ, each with amplitude pbðxÞ 5 ℱ fpðr; 0Þg: ðN (5) pbðxÞJ0 ðxrÞx dx pðr; 0Þ 5 0

This is referred to as the angular spectrum approach (Vyas and Christensen 2012; Zeng and McGough 2008). Note that a cylindrical Fourier transform, also known as a Hankel transform, is used. The function J0 ðxrÞ is the zeroth-order Bessel function of the first kind. The angular spectrum approach proceeds to calculate the propagation of each individual plane wave in the z-direction. The magnitude of the wavenumber k 5 u=c is fixed by the radial frequency of the transducer, u; and the speed of sound of the medium, c; therefore, the z-component of the wavenumber vector for each component plane wave must be given by x2z 5 k2 2x2 . If the wave were to propagate freely in an infinite fluid medium, then the pressure caused by an individual plane wave at some point z is given by px ðzÞ 5 pbðxÞexpð2ixz ðz2z0 ÞÞ. The pressure at any point caused by the superposition of all plane waves is given by

pðr; zÞ 5 5

ðN ð0N

px ðzÞJ0 ðxrÞx dx (6) pbðxÞe2ixz ðz2z0 Þ J0 ðxrÞx dx

0

The ultrasonic field of interest in this study does not propagate freely, but is confined by the components of the bioreactor. If the ultrasonic wave in the bioreactor has harmonic time dependence (radial frequency u) and the bioreactor components can be represented as a system of infinitely wide, parallel layered media (see Fig. 2b), then the extent of reflection and transmission at the interfaces between layers can be calculated using the transfer matrix method. The transfer matrix method is applied as follows. ðjÞ First, a set of scalar ð4m Þ and vectorial ðJðjÞ Þ velocity potential functions in each layer ðjÞ of the system is defined

1769

(Allard and Atalla 2009). If layer j is a fluid, then a single ðjÞ ðjÞ velocity potential function 41 (with vðjÞ 5 2V41 ) describes the ultrasonic wave. In solids, on the other hand, both a scalar velocity potential and a vectorial velocity potential are required to describe the propagation of a linear elastic wave ðjÞ ðvðjÞ 5 V41 1V3JðjÞ ; V3being the curlÞ: The vectorial ðjÞ potential JðjÞ 5 ½0; 42 ; 0 has only one non-zero component, as the other components are reduced to zero by the axial symmetry of the system. Waves in poro-elastic media are described by three velocity potentials, two of which are scalar and the third vectorial, with only a single ðjÞ non-zero component, 43 . The three velocity potentials in poro-elastic media are the eigenvectors of eqns (3) and (4) (Allard and Atalla 2009; Biot 1962). The ultrasonic field in each layer can be fully described if the propagation of the potential functions ðjÞ can be calculated. The velocity potential functions 4m are represented through the superposition of infinitely many plane waves in a similar fashion as the pressure field generated by the ultrasonic transducer: ðN  ðjÞ ðjÞ 4m ðr; zÞ 5 Am ðxÞe2ixz ðz2z0 Þ (7) 0  ixz ðz2z0 Þ 1BðjÞ ðxÞe ðxrÞx dx J 0 m The integrand in eqn (7) represents two traveling ðjÞ plane waves, each with complex amplitudes Am and ðjÞ Bm , m˛ð1Þ or ð1; 2Þ or ð1; 2; 3Þ, depending on the type of medium. The wavenumber vectors xðjÞ m depend only on the properties of layer j and the direction x of the initial component plane wave generated by the ultrasonic transducer. The coupling and boundary conditions described in the previous section provide a set of linear equations that is subsequently used to calculate the unknown comðjÞ ðjÞ plex amplitudes Am ðxÞ and Bm ðxÞ. In practice, the set of ðjÞ ðjÞ linear equations is solved to find Am ðxÞ and Bm ðxÞ for a ðjÞ ðjÞ finite set of discrete values xn . Once Am ðxn Þ and Bm ðxn Þ are known, Gaussian quadrature or the discrete Fourier transform is used to evaluate the integral in eqn (7). The resultant velocity potentials are used to calculate the velocities, pressures and stresses in each layer ðjÞ of the bioreactor at specific points ðr; zÞ. The combined method is called TM/ASA (Louw 2013). The computational cost of the TM/ASA is very low. The main drawback is the assumption that the bioreactor must be treated as a system of infinitely wide, parallel layers, as illustrated in Figure 2b. Finite volume method The TM/ASA cannot fully account for the complex geometry of the bioreactor. A more accurate solution is obtained using the FVM. FVMs are used to approximate a solution to a set of partial differential equations by

1770

Ultrasound in Medicine and Biology

integrating the equations over a small volume element (domain Uk ), assuming that all variables within the volume element are constant (for instance, pðxk Þzpk cxk ˛Uk ) (Beers 2006). The methods are very well established, and only a cursory treatment of the FVM for ultrasonic systems is provided. In the case of the sonobioreactor, the governing equations are integrated over an annular volume element with domain Umn 5 ðrm 2Dr=2; rm 1Dr=2Þ 3ðzn 2Dz=2; zn 1Dz=2Þ3ð0; 2pÞ: As an example, eqn (2) is integrated over Umn to yield

dsrr ðrm ; zn Þzðll 12GÞ dt

Volume 41, Number 6, 2015

the electrical power input, as well as the transducer output mechanical impedance Zmo , which is the ratio of the force applied at the transducer face to its velocity. The KLM transducer model (Krimholtz et al. 1970), which is based on the thickness mode vibrations of a piezoelectric slab, can be used to describe the electromechanical coupling. The KLM model is straightforward, but an in-depth explanation is beyond the scope of this article. The interested reader is referred to the original text (Krimholtz et al. 1970), as well as the excellent explanation by Bigelow (2001). Suffice to say that if the input voltage V0

ð vr ð rm 1Dr ; zn Þ ðDr12rÞ1vr ð rm 2Dr ; zn Þ ðDr22rÞ Þ 2 2 2rDr

Applying a similar approach to each of the governing equations (eqns 1–4) reduces the coupled partial differential equations to a set of ordinary differential

1ll

ð vz ð rm ; zn 1Dz Þ 2vz ð rm ; zn 2Dz ÞÞ 2 2 Dz

(8)

and the output mechanical impedance Zmo are known, then the complex velocity of the transducer face v0 can be calculated using the equations

!
  Zp 2Zmr e2ikp dp =2 2eikp dp =2 2Zp      v0 5  V0 Zei Nf Zp 1Zmf Zp 1Zmr eikp dp 2 Zp 2Zmf Zp 2Zmr e2ikp dp 

equations (Peiffer 1997), typified by eqn (8). The ordinary differential equations are subsequently reduced to an algebraic system by assuming harmonic time dependence: (vt /iu). The FVM produces a very large set of linear equations that must be solved to find the stress and velocity fields. Thus, although the FVM is suitable to the complex bioreactor geometry, it has a much higher computational cost than the TM/ASA.

Zmf 5 Zw

Zmt 5 Zp

(10)

     ! Zp 1iZmf tan kp dp 2 Zp 1iZmr tan kp dp 2  1     Zmf 1iZp tan kp dp 2 Zmr 1iZp tan kp dp 2 (11) Zei 5

KLM model of the ultrasonic transducer It is important to understand the relation between the electrical power supplied to the transducer and the mechanical power delivered to the bioreactor. The power delivered to the ultrasonic bioreactor depends on

Zmo 1iZw tanðkw dw Þ Zw 1iZmo tanðkw dw Þ

(9)

Zmt 1 1iXp 1 iuCp Nf2

(12)

According to the KLM model, all of the model parameters used in eqns (9–12) (Table 1) can be derived from the transducer properties (also given in Table 1). The transducer output mechanical impedance Zmo , on

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

1771

Table 1. KLM model parameters and transducer properties (based on a ½-in.-diameter immersion transducer, Olympus NDT, Waltham, MA, USA) KLM model parameters Active element Electrical reactance

Transducer properties

X1 5

h2 sinðkp dp Þ Zp u2

Cp 5 εAp =dp

Electrical capacitance

uZp hsinðkp dp =2Þ

Density

rp

20; 000 kg m23

Elastic stiffness

cD;p

271111i GPa

Dielectric constant

ε ε0

291

Power conversion factor

Nf 5 2

Wavenumber

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kp 5 u rp =cD;p

Electromechanical constant

h

33 V nm21

Characteristic mechanical impedance Wear plate Wavenumber

pffiffiffiffiffiffiffiffiffiffiffiffi Zp 5 Ap rp cD;p

Thickness

dp

l5 MHz =2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kw 5 u rw =cD;w

Density

rw

4000 kg m23

Characteristic mechanical impedance

pffiffiffiffiffiffiffiffiffiffiffiffiffi Zw 5 Ap rw cD;w

Elastic stiffness Thickness

cD;w dw

0:1 GPa l5 MHz =4

Backing impedance Active surface area

Zmr Ap

77 N s m21 1:27 cm2

Other

KLM 5 Krimholtz–Leedom–Matthaei (KLM) transducer model (Krimholtz et al. 1970). According to the manufacturer, the piezoelectric element and wear plate are approximately l5 MHz =2 and l5 MHz =4 thick, respectively, where l5 MHz is the ultrasonic wavelength in water at a frequency of 5 MHz:

the other hand, is not a property of the transducer but of the applied load (in this case, the bioreactor). Zmo must be known to calculate the velocity amplitude v0 of the transducer face (eqn 9). However, according to the definition of the mechanical impedance, v0 is required to calculate Zmo : Ð pdS (13) Zmo 5 S v0 S is the surface of the transducer face, and p is the pressure at surface element dS. To calculate Zmo , first note that the linearity of the governing equations ensures that if the velocity potential 4o ðxÞ with v0 5 1 is known, then the velocity potential for any other value of v0 is simply given by 4ðxÞ 5 v0 4o ðxÞ. The same is true for all other acoustic variables, such as pressure and stress. Using pðxÞ 5 v0 po ðxÞ, with po ðxÞ being the pressure calculated using v0 5 1, eqn (13) is simplified to Ð v0 po dS ð S 5 po dS (14) Zmo 5 v0 S

Coupling electrical input to ultrasonic field Applying the FVM or the TM/ASA to a specific bioreactor setup with v0 5 1 allows calculation of the pressure po at the transducer face and, subsequently, the output mechanical impedance Zmo using eqn (14). With Zmo and the applied voltage V0 as inputs, the KLM model

is used to calculate the true velocity of the transducer face, v0 (eqns 9–12). The linearity of the problem then allows the value of any acoustic variable to be determined by cðxÞ 5 v0 co ðxÞ, where co ðxÞ is calculated using v0 5 1 and cðxÞ may represent the acoustic velocity, pressure, stress, and so on. However, v0 is linearly proportional to the supplied voltage, so that cðxÞ 5 aV0 co ðxÞ, where a 5 v0 =V0 (see eqn 9). It is useful to define the normalized acoustic variables cðxÞ 5 co ðxÞ= maxjco ðxÞj as well as the gain with respect x to voltage Ac 5 maxjco ðxÞjðv0 =V0 Þ. By use of the definix tions of cðxÞ and Ac , the actual acoustic variable is calculated using cðxÞ 5 Ac cðxÞV0

(15)

In this article, the normalized acoustic variables cðxÞ and the gains Ac are reported; the combination of these two variables allows for quantitative prediction of the ultrasonic field in the sonobioreactor. The methodology is summarized in Figure 3. Experimental measurements The transducer material properties necessary for application of the KLM model are unknown and must be determined experimentally. The immersion ultrasonic transducer used in the bioreactor (5-MHz center frequency; Olympus NDT, Waltham, MA, USA) is characterized by measuring the input current and voltage and subsequently calculating the input electric impedance Zei . The experimentally measured impedance Zei is

1772

Ultrasound in Medicine and Biology

Fig. 3. Mathematical modeling methodology. Using the bioreactor parameters and assuming the transducer face velocity is fixed at v0 5 1 m=s, the scaled acoustic pressure po ðxÞ; as well as the normalized acoustic variables pðxÞ; sðxÞ; 4ðxÞ; JðxÞ; .; are calculated. The scaled pressure po ðxÞ is used to determine the transducer input mechanical impedance Zmo which, in conjunction with the Krimholtz–Leedom–Matthaei (KLM) transducer model (Krimholtz et al. 1970) and the electrical input, is used to calculate the acoustic gains Ap ; As ; A4 ; AJ and so on The normalized acoustic variables and the gains are combined to quantitatively predict the ultrasonic field in the bioreactor.

compared with the theoretically calculated impedance (using eqns 10–12), and the material properties are adjusted until the squared error is minimized. A Hewlett Packard (Palo Alto, CA, USA) 33120 A function generator was used to create a sinusoidal wave with varying frequency. The voltage drop across a resistor with known resistance was used to calculate the current i0 . Measurements were made by pre-amplifying (Panametrics 5072 PR pre-amplifier, Olympus NDT) and recording (PCSU1000 two-channel oscilloscope, Velleman, Fort Worth, TX, USA) the electrical signal. The ultrasonic transducer was kept in air, such that the output mechanical impedance is negligible and the transducer acts as a free resonator (Zmo z0). The theoretical and measured input electrical admittance Yei 5 1=Zei ˛ℂ is illustrated in Figure 4(a), and the model parameters are given in Table 1. The ratio of the

Volume 41, Number 6, 2015

Fig. 4. Measurement and modeling of ultrasonic transducer characteristics. (a) The transducer electrical admittance magnitude as predicted by the Krimholtz–Leedom–Matthaei (KLM) transducer model (Krimholtz et al. 1970) for a transducer in air. Inset: Predicted admittance for a transducer immersed in water: No wear plate resonance is seen. (b) Ratio of the velocity amplitude of the transducer face v0 ðm=sÞ to input voltage V0 ðVÞ, as a function of output mechanical impedance jZmo j and frequency. When jZmo j is small, v0 =V0 is maximized near 5 MHz, as per the transducer design. The output velocity to voltage ratio decreases when jZmo j increases.

transducer face velocity v0 to the supplied voltage as a function of frequency and output mechanical impedance is illustrated in Figure 4(b). The ultrasonic transducer is designed to be operated near 5 MHz; note the resonance effect when f 5 5 MHz (see Fig. 4a). The resonance is completely removed when the transducer is immersed in water (see inset in Fig. 4a). At an output mechanical impedance of 180 N s m21 (corresponding to the impedance created by water if no standing wave fields are

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

1773

Table 2. Sonobioreactor component properties and relevant boundary conditions Transducer Diameter Boundary conditions

12:7 mm Transducer face: velocity specified, v0 5 1 m s21 Casing: steel, approximated as rigid boundary

Tissue culture well Diameter Material properties (polystyrene) Scaffold Diameter Material properties

rf 993

Culture medium Material properties (water) Boundary conditions Hydrophone Element diameter Casing diameter Boundary conditions

34 mm (6-well plate) or 15:95 mm (24-well plate) E G rs 1,050 3,200 2,000 Kf 2,200

rs 1,163

30; 15:95 or 5 mm E G 120 46

Kb 0.11

f 0.95

G 2.0

Kf rf 993 2,200 Medium–air interface: approximated as pressure release boundary 500 mm 2:5 mm Element: impedance of water (acoustically matched) Casing: steel, approximated as rigid boundary

rs 5 Solid density (kg m-3); E 5 Young’s modulus (MPa); G 5 shear modulus (MPa); rf 5 fluid density (kg m–3); Kf 5 fluid bulk modulus (MPa); Kb 5 drained bulk modulus (MPa); f 5 porosity; G 5 tortuosity.

present), the predicted velocity amplitude is approximately 0:030 m s21 if the transducer is driven by an input voltage of 6 V (peak-to-peak) at 5 MHz. By use of the Rayleigh-Sommerfeld integral, this results in a power output of approximately 75 mW in the free field (Kinsler et al. 1999), which is typical of ultrasonic bioreactors. RESULTS Qualitative evaluation of the ultrasonic field structure The ultrasonic field is investigated using (i) scaffolds with a diameter of 5, 15.95 or 30 mm; (ii) tissue culture wells with a diameter of 15.95 mm (24-well plate) or 34 mm (6-well plate); and (ii) a hydrophone with a diameter of 2.5 mm. The simulated hydrophone is modeled on the HNR-0500 needle hydrophone (Onda, Sunnyvale, CA, USA). The ultrasonic transducer is modeled on a 12.7-mm (1/2-in.) immersion transducer by Olympus NDT. In actual experimental setups, the transducer is typically not set in an infinite baffle, as is assumed in the TM/ASA. Raising the transducer above the bottom of the bioreactor may influence the radiation pattern; the effects thereof are investigated. The results generated by the FVM and the TM/ASA are qualitatively compared to determine whether the TM/ ASA captures the dominant features of the ultrasonic field. The normalized pressure field pðxÞ in the bioreactor is reported in all figures (in solid and porous layers, the normalized total normal stress sðxÞ is shown instead of the pressure); this can readily be converted to the actual pressure if the gain Ap and the supplied voltages are known (eqn 15). Table 2 summarizes the material properties and geometries of the bioreactor components. The relevant

boundary conditions used in the numerical simulations are also given in Table 2. A non-reflecting boundary condition is always applied at the radial boundary when using the FVM. Ultrasonic field in a 6-well-plate, 30-mm-diameter scaffold Figure 5 compares the results generated by (a) the TM/ASA, (b) the FVM with a transducer set in an infinite, rigid baffle and (c) the FVM with the transducer raised above the bioreactor bottom. In each case, the diameter of the scaffold (30 mm) is greater than that of the ultrasonic transducer (12.7 mm). The formation of a standing wave field is immediately obvious, with regularly spaced nodes and antinodes in the water layers and the well bottom. The standing wave field in the water qualitatively matches experimental measurements reported in the literature (Leskinen and Hynynen 2012). The model clearly predicts a similar standing wave field in the porous scaffold, the presence of which will affect tissue engineering experiments. The ultrasonic beam is highly directive. The TM/ ASA predicts that all of the acoustic pressure is focused in a single beam directly ahead of the transducer; no secondary lobes are formed. The FVM, on the other hand, predicts a secondary lobe with a maximum pressure along a 4 line drawn from the edge of the transducer. When the transducer is raised (i.e., not set within an infinite baffle), the maximum pressure of the secondary lobe occurs at an angle of approximately 7 . The amount of acoustic power contained in the side lobes is negligible compared with that in the primary beam. Once again, these predictions are in line with reported experimental measurements (Leskinen and Hynynen 2012).

1774

Ultrasound in Medicine and Biology

Volume 41, Number 6, 2015

Fig. 5. Normalized ultrasonic pressure field in a bioreactor with a large-diameter scaffold. (a) Transfer matrix/angular spectrum approach (TM/ASA)-predicted pressure field. A single beam is visible, with a clear standing wave field. (b) Finite volume method (FVM)-predicted pressure field for a transducer set in an infinite, rigid baffle. Side lobes are formed, the first at an angle of 4 , but the pressure in the side lobes is negligible compared with that in the primary beam (c) FVM-predicted pressure field for a raised transducer. Raising the transducer increases the angle of the first side lobe to 7 . The ultrasonic field exhibits high pressure bands in the scaffold in all three cases. The pressure is maximized on the axis r 5 0.

The ultrasonic field in the water layer directly adjacent to the transducer exhibits an even pressure distribution in the radial direction when r,a (with a being the radius of the transducer), because of the high beam directivity. However, the mathematical model predicts highand low-pressure bands within the porous scaffold, with a clear high-pressure band along the central axis of the transducer. The formation of high-pressure bands can be attributed to interference by the various wave modes present in the scaffold: two longitudinal waves and a single shear wave. The pressure striations will result in unequal mechanical stimulation of cells seeded throughout the scaffold. Effects of decreasing scaffold and well diameters Figure 6(a) illustrates the FVM-predicted ultrasonic field when a 5-mm-diameter scaffold is placed in a 34-mm-diameter well. Qualitatively, the structure of the ultrasonic field becomes much more complex. Low-pressure zones are formed around the edges of the scaffold, and the pressure field within the scaffold is concentrated in the middle, resulting in uneven mechan-

ical stimulation of cells. This effect is repeated when the scaffold is placed in a 15.95-mm-diameter well (Fig. 6b). If, however, a 15.95-mm-diameter scaffold is placed in the 15.95-mm-diameter well, the pressure field reverts back to an even, albeit slightly banded, structure. The pressure distribution in the 15.95-mm-diameter scaffold is comparable to that of the 30-mm-diameter scaffold, as can be seen by comparing Figures 5 and 6(c). In this case, all acoustic elements have a diameter at least 25% greater than that of the ultrasonic beam. Thus, radial reflections from the well wall are negligible. Leskinen and Hynynen (2012) experimentally measured surface waves with a velocity amplitude of 0.030 m/s originating from the well wall when operating a transducer at 1.035 MHz. The velocity amplitude decreased to 0.005 m/s when the transducer was operated at 3.35 MHz. The decrease is explained by noting that an increase in frequency would correspond to higher beam directivity and, thus, reduced wall effects. The effects of the well wall will therefore become negligible at even higher frequencies, as is the case in the simulations above,

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

1775

Fig. 6. Normalized ultrasonic pressure field in the bioreactor when a scaffold is used with a diameter smaller than that of the transducer. (a) Pressure field given a scaffold with diameter 5 mm, set in a well with diameter 34 mm. The pressure becomes concentrated toward the center of the scaffold, with low pressure zones being generated around the edges. (b) Pressure field for the same scaffold set in a well with diameter 15:95 mm. The narrower well has very little effect on the ultrasonic field. (c) If the 15.95-mm-diameter well is filled with a scaffold, a more evenly distributed pressure field is observed throughout the scaffold, once again comparable to transfer matrix/angular spectrum approach (TM/ASA) results. (d) A hydrophone with a diameter of 2.5 mm interferes with the pressure distribution, creating further concentration of pressure on the central axis.

which predict the ultrasonic field generated by a transducer operating at 5 MHz. It is important to compare the diameter of the bioreactor components with the diameter of the ultrasonic beam, not the diameter of the transducer, when considering radial reflections and wall effects. Hydrophone effects The hydrophone investigated here is obtrusive, as the diameter of the hydrophone casing is 2.5 mm (even though the active element is only 0.5 mm). Figure 6(d) illustrates that the hydrophone further concentrates the pressure near the central axis. It is therefore important to account for the hydrophone geometry when experimentally measuring the pressure distribution. The effect of the hydrophone on the ultrasonic field will become negligible if a hydrophone is used with dimensions substantially smaller than the acoustic wavelength. In the current setup, lz300 mm, thus necessitating the use of fiberoptics-based hydrophones (e.g., HFO series, Onda) or similar.

Quantitative analysis of different bioreactor setups Besides the qualitative structure of the ultrasonic field (described by the normalized pressure pðxÞ and normalized stress sðxÞ as discussed above, it is important to identify quantitative indicators of the mechanical stimulation delivered to cells in a bioreactor. The primary factors affecting the applied stimulation were (i) resonance effects and (ii) nodal positions. The extent of these effects was investigated in a variety of different bioreactor setups, as listed in Table 3. Resonance effects The bottom-up configuration (Table 3, reactor setups 1 and 3) is sensitive to changes in the liquid level. The liquid level sensitivity has been experimentally observed, causing deviations of up to 700% (Leskinen and Hynynen 2012). Figure 7(a) illustrates the sensitivity of the pressure to liquid level. If the transducer face velocity v0 is fixed, the maximum pressure gain with respect to the

1776

Ultrasound in Medicine and Biology

Volume 41, Number 6, 2015

Table 3. Various sonobioreactor setups Setup

Layer material (thickness, mm)*

1. Bottom-up with scaffold

jj Water (10)/polystyrene (1.52)/scaffold (5)/water (3)/air (0.1) jj Water (10)/polystyrene (1.52)/water (8)/air (0.1) Air (0.1)/polystyrene (1.52)/scaffold (5)/water (3) jj Air (0.1)/polystyrene (1.52)/water (3) jj

2. Bottom-up without scaffold 3. Top-down w scaffold 4. Top-down w/o scaffold

* jj marks transducer position.

transducer face velocity (i.e., maxjpðxÞj=v0 ) varies bex tween 2000 and 4500 kPa/m s21. However, the pressure variation in the bioreactor influences the output mechanical impedance of the transducer and, subsequently, the ratio of the transducer face velocity to applied voltage, v0 =V0 (eqns 9–12). If the supplied voltage were fixed instead of the transducer face velocity, the maximum pressure gain as a function of liquid level would vary between 20 and 180 kPa/V. The resonance peaks are also much sharper, as illustrated in Figure 7(a). The difference in functional form between the pressure gain with v0 fixed and the pressure gain with V0 fixed is due to the inclusion of the changing transducer output mechanical impedance. Figure 7(a) highlights the necessity to couple the ultrasonic field model to an appropriate transducer model (such as the KLM model). A correlation exists between the maximum pressure amplitude and the input electrical impedance, as can be seen in Figure 7(a, b). The input electrical impedance, Zei , which is readily measured, can be used to infer the pressure amplitude without the use of a hydrophone. If the impedance begins to change, the input voltage, V0 ; can be lowered to limit pressure peaks. The major difference between bottom-up and top-down sonication is that the distance between the transducer and the air/water interface can vary with evaporation or immersion of objects during bottom-up sonication, but the distance between the transducer and the reactor bottom/air interface remains constant for the top-down configuration. The top-down configuration is subsequently not vulnerable to resonances if the liquid level changes slightly, but will be sensitive to resonance when the applied frequency is adjusted. Nodal positions In all setups, the ultrasonic wave is almost completely reflected by the air interface because of the acoustic impedance mismatch, resulting in a standing wave field with clear pressure nodes and antinodes. Previous work has indicated that cells respond differently to different forms of mechanical stimulation, that is, dila-

Fig. 7. Pressure gain as a function of liquid level. (a) The pressure gain varies with liquid level when the transducer face velocity, v0 ; is fixed. However, the pressure in the bioreactor is dependent on the output mechanical impedance, Zmo , which in turn depends on the liquid level in the tissue culture well. The maximum pressure gain as a function of liquid level is very different when the input voltage is fixed, as opposed to the transducer face velocity. (b) The input electrical impedance is also dependent on the liquid level in the bioreactor. If the input electrical impedance is used as a measured variable, the pressure amplitude can be modulated by adjusting the input voltage.

tion or shear. Furthermore, cells in standing ultrasonic fields will experience either shear or dilation depending on their proximity to pressure nodes and antinodes (Louw et al. 2013). During bottom-up sonication, all pressure nodes throughout the system are spaced a distance nl=2 from the air–water interface, where n˛Z and l is the wavelength. Therefore, the positions of pressure nodes depend on the liquid level above the scaffold, as well as the frequency at which the transducer is operated. During top-down sonication, the pressure node and antinode positions likewise depend on the frequency of the applied ultrasound, but not the liquid level. Figure 8

Evaluation of acoustic field in ultrasonic bioreactor d T. M. LOUW et al.

1777

configuration is ideal for investigating the effect of a very specific type of mechanical stimulation on cells. DISCUSSION

Fig. 8. Effect of frequency on the positions of pressure nodes in bioreactor setup 4. When no scaffold is present, cells will be seeded on the bottom of the tissue culture well, indicated above by a red line. The positions of pressure nodes and antinodes depend on the frequency of the applied ultrasound. When f 5 3 MHz, cells will be located at a pressure node (minimum acoustic pressure), whereas they will be located at a pressure antinode (maximum acoustic pressure) if f 5 5 MHz.

illustrates the normalized pressure distribution near the well bottom when bioreactor setup 4 is used (top down, no scaffold). At 3 MHz, the cell layer lies at a pressure node and cells would experience shear stress. At 5 MHz, the cell layer lies at or near pressure anti-nodes, and the mechanical stimulation would be mostly dilational strain. At 4 MHz, the cell layer lies between a node and an antinode; thus, a mixture of shear and dilational strain can be expected. Therefore, the top-down

The TM/ASA method is a semi-analytical technique that can be used to analyze different bioreactor layouts. The analysis above indicates that the TM/ASA produces results comparable to those of the FVM whenever each acoustic component has a diameter at least 25% larger than that of the ultrasonic beam, in which case negligible radial reflection will occur. In all other cases, the FVM should be used to predict the ultrasonic field. The main advantage of the TM/ASA is the low associated computational cost, thus allowing investigation of the effects of continuously varying parameters of interest (such as the liquid level) (Fig. 7). The model is especially useful for investigating the ultrasonic field in regions where measurements are impossible, such as in the porous scaffold. Within the scaffold, the model predicts the formation of vertical high-pressure bands in addition to the expected pressure variations associated with standing wave nodes and antinodes. The banded structure of the ultrasonic field will undoubtedly influence the mechanical stimulation experienced by cells seeded on the scaffold, thereby compromising the mechanical integrity of any synthetic tissue. The ultrasonic field becomes even more irregular when the scaffold diameter is smaller than that of the ultrasonic beam. It is recommended that the diameter of the scaffolds in tissue engineering experiments be at least 25% larger than the diameter of the transducer. Besides the ultrasonic field structure, the main sources of experimental error are variations in liquid levels (Fig. 7) and nodal positions (Fig. 8). The results highlight the necessity of coupling any acoustic model to an appropriate transducer model to recreate resonance effects. Errors caused by resonance effects may be avoided by using the TM/ASA–KLM coupled model to design experiments correctly while concurrently using voltage and current measurements at the transducer electrical input to infer the pressure amplitude in the bioreactor. Although the ultrasonic field is sensitive to external disturbances, repeatable experiments can be designed by understanding the primary factors affecting the pressure amplitude, nodal positions and beam shape; mathematical models provide the necessary insights. CONCLUSIONS Experimental studies have unveiled the complexity of ultrasonic bioreactors and have indicated the need for in-depth analyses. The TM/ASA partly addresses this need: It is a useful technique for evaluating the effects of various factors (such as the liquid level) because of its

1778

Ultrasound in Medicine and Biology

low computational cost. However, it is limited to specific reactor geometries where radial reflections can be neglected; in all other cases, the FVM is more appropriate. The mathematical model further allows investigation of the ultrasonic field within the porous scaffold. Additional work is required to experimentally validate the model and investigate parameter sensitivity. However, the qualitative insights gained will provide much needed guidance to experimentalists in ultrasound-enhanced tissue engineering. Acknowledgments—We acknowledge Kelly Hudek for conducting experimental work on transducer characterization.—This work was supported, in part, by American Recovery and Reinvestment Act of 2009 Research Grant 1 R21 RR024437-01 A1 from the Department of Health and Human Services.

SUPPLEMENTARY DATA Supplementary data related to this article can be found online at http://dx.doi.org/10.1016/j.ultrasmedbio.2014.12.015.

REFERENCES Allard J, Atalla N. Propagation of sound in porous media. 2nd ed. New York: Wiley; 2009. Beers KJ. Boundary value problems. In: Numerical methods for chemical engineering: Applications in MATLAB. Cambridge, MA: Cambridge University Press; 2006. p. 258–316. Bigelow TA. Experimental evaluation of non-linear indices for ultrasound transducer characterizations. Thesis, University of Illinois at Urbana–Champaign, 2001. Biot MA. Generalized theory of acoustic propagation in porous dissipative media. J Acoust Soc Am 1962;34:1254. Dalecki D. Mechanical bioeffects of ultrasound. Annu Rev Biomed Eng 2004;6:229–248. Fowlkes JB. American Institute of Ultrasound in Medicine consensus report on potential bioeffects of diagnostic ultrasound: Executive summary. J Ultrasound Med 2008;27:503–515. Hasanova GI, Noriega SE, Mamedov TG, Guha Thakurta S, Turner JA, Subramanian A. The effect of ultrasound stimulation on the gene and protein expression of chondrocytes seeded in chitosan scaffolds. J Tissue Eng Regen Med 2011;5:815–822. Hensel K, Mienkina MP, Schmitz G. Analysis of ultrasound fields in cell culture wells for in vitro ultrasound therapy experiments. Ultrasound Med Biol 2011;37:2105–2115. Hsu HC, Fong YC, Chang CS, Hsu CJ, Hsu SF, Lin JG, Fu WM, Yang RS, Tang CH. Ultrasound induces cyclooxygenase-2 expression through integrin, integrin-linked kinase, Akt, NF-kappaB and p300 pathway in human chondrocytes. Cell Signal 2007;19:2317–2328. Hsu S, Huang T, Chuang S, Tsai IJ, Chen DC. Ultrasound preexposure improves endothelial cell binding and retention on biomaterial surfaces. J Biomed Mater Res B Appl Biomater 2006a;76: 85–92. Hsu S, Kuo CC, Whu SW, Lin CH, Tsai CL. The effect of ultrasound stimulation versus bioreactors on neocartilage formation in tissue engineering scaffolds seeded with human chondrocytes in vitro. Biomol Eng 2006b;23:259–264.

Volume 41, Number 6, 2015 Karande P, Jain A, Ergun K, Kispersky V, Mitragotri S. Design principles of chemical penetration enhancers for transdermal drug delivery. Proc Natl Acad Sci USA 2005;102:4688–4693. Kinoshita M, Hynynen K. Key factors that affect sonoporation efficiency in in vitro settings: The importance of standing wave in sonoporation. Biochem Biophys Res Commun 2007;359:860–865. Kinsler LE, Frey AR, Coppens AB, Sanders JV. Fundamentals of acoustics. 4th ed. New York: Wiley; 1999. Krimholtz R, Leedom DA, Matthaei GL. New equivalent circuits for elementary piezoelectric transducers. Electron Lett 1970;6:398. Leskinen JJ, Hynynen K. Study of factors affecting the magnitude and nature of ultrasound exposure with in vitro set-ups. Ultrasound Med Biol 2012;38:777–794. Louw TM. Mathematical modeling of ultrasound in tissue engineering: From bioreactors to the cellular scale. Thesis, University of Nebraska-Lincoln, 2013. Louw TM, Budhiraja G, Viljoen HJ, Subramanian A. Mechanotransduction of ultrasound is frequency dependent below the cavitation threshold. Ultrasound Med Biol 2013;39:1303–1319. Mitragotri S. Healing sound: The use of ultrasound in drug delivery and other therapeutic applications. Nat Rev Drug Discov 2005;4: 255–260. Naito K, Watari T, Muta T, Furuhata A, Iwase H, Igarashi M, Kurosawa H, Nagaoka I, Kaneko K. Low-intensity pulsed ultrasound (LIPUS) increases the articular cartilage type II collagen in a rat osteoarthritis model. J Orthop Res 2010;28:361–369. Nelson TR, Fowlkes JB, Abramowicz JS, Church CC. Ultrasound biosafety considerations for the practicing sonographer and sonologist. J Ultrasound Med 2009;28:139–150. Nishikori T, Ochi M, Uchio Y, Maniwa S, Kataoka H, Kawasaki K, Katsube K, Kuriwaka M. Effects of low-intensity pulsed ultrasound on proliferation and chondroitin sulfate synthesis of cultured chondrocytes embedded in Atelocollagen gel. J Biomed Mater Res 2002; 59:201–206. Noriega S, Hasanova G, Subramanian A. The effect of ultrasound stimulation on the cytoskeletal organization of chondrocytes seeded in three-dimensional matrices. Cells Tissues Organs 2013;197:14–26. Parvizi J, Wu CC, Lewallen DG, Greenleaf JF, Bolander ME. Low-intensity ultrasound stimulates proteoglycan synthesis in rat chondrocytes by increasing aggrecan gene expression. J Orthop Res 1999; 17:488–494. Peiffer A. The acoustic finite integration technique for waves of cylindrical symmetry (CAFIT). J Acoust Soc Am 1997;102:697. Rubin C, Bolander M, Ryaby JP, Hadjiargyrou M. The use of lowintensity ultrasound to accelerate the healing of fractures. J Bone Joint Surg Am 2001;83A:259–270. Saenger EH, Gold N, Shapiro SA. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion 2000; 31:77–92. Virieux J. P-SV wave propagation in heterogeneous media: Velocity stress finite difference method. Geophysics 1986;51:889–901. Vyas U, Christensen D. Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method. IEEE Trans Ultrason Ferroelectr Freq Control 2012;59:1093–1100. Whitney N, Lamb A. Integrin-mediated mechanotransduction pathway of low-intensity continuous ultrasound in human chondrocytes. Ultrasound Med Biol 2012;38:1734–1743. Zeng X, McGough RJ. Evaluation of the angular spectrum approach for simulations of near-field pressures. J Acoust Soc Am 2008;123: 68–76. Zhang ZJ, Huckle J, Francomano CA, Spencer RGS. The effects of pulsed low-intensity ultrasound on chondrocyte viability, proliferation, gene expression and matrix production. Ultrasound Med Biol 2003;29:1645–1651.