Computational and Mathematical Methods in Medicine Computational Methods and Models in Circulatory and Reproductive Systems

The circulatory and reproductive systems can be considered as tubular transport systems of blood cells and gametes, respectively. Both systems involve complex flows and fluid-structure interactions of large displacements and large deformations. In addition, multiple scales and multiple processes need to be considered. Because of their importance in health and fundamentals, great effort has been put in developing numerical methods for these systems and in understanding the underlying flow physics. In spite of the many achievements, it is still a challenge to model these systems accurately using numerical methods. Moreover, there still remains much concerted development towards grasping the complicated flow physics within such systems by using advanced numerical methods and mathematical models. Therefore, the goal of this special issue was to collect together recent contributions on the development of numerical methods for complex flows and fluid-structure interaction in circulatory and reproductive systems and/or on addressing the complicated flow physics and heat and mass transfer in these systems for fundamental understanding as well as engineering applications. The papers can be divided into three groups as detailed below.


Introduction
The circulatory and reproductive systems can be considered as tubular transport systems of blood cells and gametes, respectively. Both systems involve complex flows and fluidstructure interactions of large displacements and large deformations. In addition, multiple scales and multiple processes need to be considered. Because of their importance in health and fundamentals, great effort has been put in developing numerical methods for these systems and in understanding the underlying flow physics. In spite of the many achievements, it is still a challenge to model these systems accurately using numerical methods. Moreover, there still remains much concerted development towards grasping the complicated flow physics within such systems by using advanced numerical methods and mathematical models. Therefore, the goal of this special issue was to collect together recent contributions on the development of numerical methods for complex flows and fluid-structure interaction in circulatory and reproductive systems and/or on addressing the complicated flow physics and heat and mass transfer in these systems for fundamental understanding as well as engineering applications. The papers can be divided into three groups as detailed below.

Patient-Specific Geometry Simulations
The work by J. S. Byun, S.-Y. Choi, and T. Seo attempts to investigate the hemodynamic phenomena in the cerebral arteries before and after surgery of the aneurysm with patient-specific pre-and postsurgery cerebral arterial geometries by using three-dimensional computational fluid dynamics models. In the models, the flow is approximated by an incompressible and Newtonian fluid in laminar flow regime, and the vascular wall is assumed to be rigid. The flow patterns, the inflow jet streams, and the wall average shear stress have been considered by using different patient-specific geometries. Several interesting findings are reported, among which the average shear stress distribution associated with aneurysm rupture is attractive and could provide useful insight in predicting the risk of aneurysm rupture.
X. Liu et al. present a numerical study of flow characteristics in the patient-specific upper airway obstructed by the pharyngeal collapse. In this study, the patient-specific geometries are used, the flow dynamics is solved by using computational fluid dynamics, and the turbulence is modeled by the large eddy simulation. The numerical simulations are validated with their measurement by a laser Doppler anemometry. Simulations are conducted by varying flux rate under both continuous inspiration and expiration with focus on discussing the velocity fields and static pressure fields and their correlation with the upper airway statuses (e.g., narrowing of pharynx). Such information could provide guidance for the treatment of obstructed respiratory disease. In addition, the numerical procedures reported can be directly applied to study flows in narrowed artery, arteriovenous graft, and arteriovenous fistula.

Cell/Capsule Dynamics
J.-T. Ma, Y.-Q. Xu, and X.-Y. Tang report a size-dependent cell sorting scheme based on a controllable asymmetric pinched flow by using an immersed boundary-lattice Boltzmann method. The scheme employs a device which consists of 2 upstream branches, 1 transitional channel, and 4 downstream branches. To study the cell sorting efficiency, systematic simulations are conducted by varying the inlet flow ratio, the cell size, and the outlet flow ratio. The device is approved to be effective by the finding that cells of different diameters can be successfully sorted into different downstream branches. This work provides a reference for the design of microfluidics for cell/particle sorting.
The work by F.-B. Tian presents an immersed boundarylattice Boltzmann method for capsule fluid-structure interaction involving non-Newtonian fluid (e.g., power-law fluid). In this method, the capsule dynamics and the fluid dynamics are coupled by using the immersed boundary method; the incompressible viscous power-law fluid dynamics is obtained by solving the lattice Boltzmann equation where a shear rate-dependant relaxation time is used to achieve the non-Newtonian rheology. The non-Newtonian flow solver is validated by considering a power-law flow in a channel and then applied to study the deformation of a capsule in a power-law shear flow by varying the Reynolds number, dimensionless shear rate, and power-law index. Several interesting results that are distinctive in the power-law flow are reported. This work provides an alternative for fluid-structure interactions involving non-Newtonian fluid. Further effort can be made to explore applications of non-Newtonian effects in cell/particle sorting.

Molecular Dynamics Simulation
Y. Ge et al. present a study on the effects of solution concentration on ion distribution in a nanopore-based device inspired from red blood cells by using a molecular dynamics method. It is found that the density peaks for both the counterion and coion near the charged wall increase at different speeds with the increase of the solution concentration if the screening effects appeared, and consequently the potential near the charged wall of the nanopore switches from a negative value to a positive one during the simulation. This work provides insights in controlling the ion permeability and improving the cell transfection as well as the design of nanofluidic devices. In addition, it inspires future work on multiscale methods including the method in this work and those methods in the studies by F.-B. Tian and J.-T. Ma, Y.-Q. Xu, and X.-Y. Tang to tackle multiscale cell dynamics.
The wide biomedical and mechanical engineering science community will be interested in this special issue. Both researchers and practitioners will find results presented in this issue helpful in many applications.

Introduction
Flow induced deformation of a capsule consisting of a membrane enclosing an internal medium such as a gel or a liquid is an important problem in fundamental research as well as bioengineering applications. For example, a capsule in shear flow is a fundamental process that is related to erythrocytes (or red blood cells), leukocytes (or white blood cells), and platelets in blood flow [1][2][3][4][5][6]. Deformation is essential for red blood cells to perform their physiological functions in the circulation of capillary blood vessels and thus affects the rheology of the blood [6][7][8]. The deformations of white blood cells and red blood cells can, respectively, affect the immune response and the oxygen load release [9,10]. The synthetic microcapsules with polymerized interfaces are designed for drug delivery, cosmetic production, and other technical usages [11,12]. Therefore, great effort has been made to study this problem (e.g., [1,4,6,8,10,[12][13][14]).
Both experimental and numerical methods have been conducted to observe capsule behaviors and the relevant underneath fluid-structure interaction physics. Schmid-Schönbein and Wells [15] and Goldsmith [16] observed that red blood cells tumble like rigid particles at low shear rates while they deform to a steady configuration and direction after which the membrane rotates around the internal liquid (tank-treading movement) at high shear rates. Later, Goldsmith and Marlow [17] and Keller and Skalak [18] found that the viscosity ratio between the liquids inside and outside the cell may also affect the type of behaviors. A higher viscosity inside would cause unsteady tumbling-rotating motion, while a smaller viscosity inside would lead to the tank-treading movement with a stationary shape. These phenomena were captured by Xu et al. [14]. More recently, Dupire et al. [19] reported rolling motion in addition to other behaviors. A hysteresis cycle and two transient dynamics driven by the shear rate (i.e., an intermittent regime during the "tanktreading-to-flipping" transition and a Frisbee-like "spinning" regime during the "rolling-to-tank-treading" transition) were highlighted.
There are several numerical methods that have been used to study capsule dynamics. Examples are the boundary element method (e.g., [20]), arbitrary Lagrangian-Euler method (e.g., [21][22][23]), immersed finite element method (e.g., [24]), and immersed boundary method (IBM) (e.g., [ [12][13][14][25][26][27][28][29][30][31][32][33][34]). Specifically, Zhou and Pozrikidis [20] studied the transient and large deformation of capsules with position-dependent membrane tension. Choi and Kim [21] simulated the motion of red blood cells freely suspended in shear flow to investigate the nature of pairwise interception of red blood cells using a fluid-particle interaction method based on the arbitrary Lagrangian-Eulerian method. Villone et al. [22,23] studied the effect of the non-Newtonian fluid on flexible particle deformation and migration in shear and channel flows by using the arbitrary Lagrangian-Eulerian method. The Navier-Stokes equations and cell-cell interaction were coupled in the framework of the immersed finite element method and mesh-free method by Y. Liu and W. K. Liu [24] to model complex blood flows with deformable red blood cells within micro and capillary vessels in three dimensions. The transient deformation of a liquid-filled elastic capsule in simple shear flow was studied by Sui et al. [1,4,35,36]. The fluid inertia on the dynamics of deformable particles has been studied by Krüger et al. [32] and Kaoui and Harting [34]. More recently, optical force based separation of particles/capsules was simulated by Chang et al. [37][38][39]. Still, as far as known to us, the existing numerical simulations seldom consider the non-Newtonian rheology effects on the capsule behaviors, while blood and most fluids involved in biomedical engineering are non-Newtonian fluids [6,8,40,41].
Following the work by Sui et al. [1] and Xu et al. [14], we develop an immersed boundary-lattice Boltzmann method (IB-LBM) to study the non-Newtonian effects on the deformation of a capsule in a shear flow. As a typical rheology, the power-law fluid is used. In this method, the capsule dynamics and the fluid dynamics are coupled by using the IBM, and the incompressible viscous power-law fluid motion is acquired by solving the lattice Boltzmann equation (LBE).
The rest of this paper is organized as follows. Section 2 briefly introduces the governing equations of the fluid and solid structures and describes the numerical approach. Section 3 presents the numerical results. Final conclusions are given in Section 4.

Physical Model and Mathematical Formulation.
In this work, a two-dimensional liquid capsule enclosed by an elastic membrane and immersed in an incompressible non-Newtonian fluid is considered, as illustrated in Figure 1 where is the arch length coordinate, n denotes the surface normal that points into the outer fluid, t denotes the tangent unit vector that points to the increasing arc length, and 0 is the velocity applied at both top and bottom walls to form a simple shear flow. The incompressible non-Newtonian fluid dynamics is achieved by using LBM [42,43]. Great effort has been made in using LBM to solve the complex flows (see several reviews [42][43][44] for the effort). Many publications have presented the details of LBM; thus we just provide a brief description in this paper and discuss the extension for non-Newtonian fluids. The details of LBM and its applications are referred to the references provided. Using the IB-LBM, the lattice Boltzmann equation (LBE) that governs the viscous flow dynamics and incorporates the traction jump across the interface due to the elastic membrane is written as [1,14,42,43,45,46] where (x, ) is the distribution function for particles with velocity e at position x and time , Δ is the size of the time step, eq (x, t) is the equilibrium distribution function, LB represents the dimensionless relaxation time, is the term representing the body force effect on the distribution function, are the weights, u = ( , V) is the velocity of the fluid, is the speed of sound defined by = Δ / √ 3Δ with Δ being grid spacing, f is the body force acting on the fluid, ΔF( , ) is the Lagrangian force density on the fluid by the elastic boundary, X is the position vector on the membrane, and (x − X(s, )) is Dirac's delta function.
In the present work, a two-dimensional nine-speed (D2Q9) model is used, as shown in Figure 2. In this model, the nine possible particle velocities are given by The values of e ensure that, within one time step, a particle moves to one of the eight neighboring nodes as shown in Figure 2 or stays at its current location. The weights, , are given by 0 = 4/9 and = 1/9 for = 1 to 4 and = 1/36  Figure 2: Nine base vectors representing 9 possible velocity directions in the D2Q9 lattice model. for = 5 to 8. In addition, the relaxation time is related to the kinematic viscosity in the Navier-Stokes equations in terms of where ] = / with being the dynamic viscosity of the ambient fluid and being the fluid density. When the particle density distributions are known, the fluid density, velocity, and pressure are then computed from = ∑ , u = ∑ e + 0.5fΔ , Theoretically the LBM introduced above simulates the compressible viscous flow instead of incompressible viscous one, because the spatial density variation is not zero in LBM simulations. In the applications, the Mach number (Ma = 0 / ) should be low (e.g., Ma ≤ 0.3) so that the incompressible viscous flow can be correctly simulated. The deduction process from LBE to the incompressible viscous flow governing equations can be found in [47]. The dynamics viscosity is a constant for a Newtonian fluid, while it is dependent on the local shear rate for a non-Newtonian fluid. Without loss of generality, the power-law fluid is taken as a representation of non-Newtonian fluids in the present paper. The rheological equation of state for powerlaw fluids is defined by [48] =̇− 1 , = max (√2 ,̇) , where is the power-law consistency index, is the powerlaw fluid behavior index,̇is the shear rate, anḋis the minimum shear rate that is applied to avoid the numerical singularity caused by the zero shear rate. The power-law fluids of < 1, > 1 and = 1 are, respectively, the shear-thinning, shear-thickening, and Newtonian fluids. In (9), the Einstein summation convention is applied. In LBM implementation, can be either calculated macroscopically by using the finite difference method or locally in mesoscopic scale by using eq and f [49]. To achieve the non-Newtonian rheology, a shear rate-dependant relaxation time is used which can be obtained by applying the effective viscosity determined by (8) in (6).
Because of the deformation, the membrane develops a transverse shear tension and a bending moment . In addition, due to the stretching motion, a tension, , is induced. Consider the force balance of membranes; we acquire Please note that ΔF is the Lagrangian force on the fluid exerted by the elastic body boundary and is opposite to the fluid force on the boundary. To evaluate and for the thin membrane, we use Hooke's law which is a relatively simpler constitutive law for modeling small deformation of capsules. Hooke's law states that the tension and the bending moment are linearly related to the stretch and the curvature, respectively. It can be written in the form where is the bending coefficient, is the stretching coefficient, 0 is the initial arch length, is the curvature of membrane, and 0 is the curvature in the resting configuration. If is large so that the stretching deformation is small, Hooke's law works well. Following the work by Sui et al. [1], the capsule membrane is assumed to be infinitely thin so that the bending effect is neglected; that is, = 0. Actually, the effect of is similar to that of when is small compared to [1,35]. If is large, the capsule may undergo flipping motion [35].
The velocity of a point on the capsule is interpolated from the flow field, and the position of the capsule is updated explicitly; that is, where U( , ) is the velocity of the capsule.
In this work, we choose the flow shear rate (e.g., 0 /ℎ), density, and the radius of the capsule to nondimension the governing equations and obtain two dimensionless parameters: the Reynolds number Re and dimensionless shear rate , which are defined by where is the radius of the undeformed capsule. measures the ratio of shear force to the elastic force. For applications where inertia force is important, we can also use Re ⋅ to nondimension the elastic property, which measures the ratio of fluid inertial forces to stretching elastic forces. Please note that the two-dimensional model is used in this work, while red blood cell deformation is a three-dimensional problem; however, the results obtained in this research should show some features common with the three-dimensional simulations, as demonstrated in [1].

Numerical
Method. Similar to [1], the capsule is discretized by nodal points which are initially distributed with equal distances. The position of the th node at time level is denoted by X . To compute the stretching force at th node, a finite difference scheme is used; that is, where Δ is the Lagrangian grid spacing on the membrane and the tension and tangent vector, t = X/ , at the segment center, + 1/2, are both computed using a secondorder central difference scheme. The time integration of (14) is calculated according to In the IBM, a smooth approximation [50] of Dirac's delta function, ℎ , is used, In the present simulations, Δ = Δ = Δ (in lattice units) is used. Now, the computational algorithm can be summarized as follows. Given all values at time step , the values at time step + 1 can be undated by the following: (1) Calculate the Lagrangian force density ΔF +1 from X by using (11)- (12).
In the present work, the above-mentioned computational simulation algorithm is implemented in the Fortran 90 programming language.

Validation.
The IB-LBM in this work has been validated and verified in our previous studies (see, e.g., [14,46]) and has been used to study filament(s) flapping in viscous fluids [51][52][53], sperm swimming, and cell/particle flows [10,54]. In the present work, we focus on the validation of non-Newtonian flow by considering a power-law flow in a straight channel which is one of the benchmark problems to validate an in-house computational fluid dynamics solver. As in our previous work [41], we consider a two-dimensional steady laminar developing flow of power-law fluid with a uniform incoming velocity ∞ through a rectangular channel of height ℎ and length , as shown in Figure 3. The physically realistic initial and boundary conditions are given as The computations are performed with the dimensionless domain size ( /ℎ × 1) of 40 × 1 discretized by 2001 × 51   Figure 4 with the corresponding analytical solution for fully developed velocity profile [41,48] for power-law fluid flow in a channel which is given as From Figure 4, it is found that the present numerical results show a good agreement with the analytical solutions for various values of power-law index, giving us confidence in the reliability and accuracy of the present numerical solution procedure. It is noted from Figure 4 that the shear layer is thinned for < 1 and thickened for > 1 compared to the Newtonian fluid case ( = 1).

Numerical Results
We first consider the power-law index effect on the deformation of a cylindrical capsule in a shear flow. The Reynolds number is 0.05, which is in the range of normal physiological conditions. The dimensionless shear rate is 0.04. The computational domain ranges from 0 to 20 in bothaxis and -axis. The capsule is at the center of the domain, and its membrane is equally discretized into 80 Lagrangian nodes. The grid resolution is Δ = Δ = Δ = 0.1 .
The characteristic velocity is set as 0 = 0.05 so that the dimensionless relaxation time is 0.5 < LB < 3.0. Such setup is consistent with that used in [1]. To study the power-law index effect, is set in the range of 0.2 < < 1.8, covering the shearthinning, Newtonian, and shear-thickening fluids. In order to quantify the deformation of the capsule, the Taylor shape parameter is introduced [1], where and B are, respectively, the length and width of a cross-section of the cylindrical capsule. Figure 5 shows the deformation of the flexible capsule in a shear flow for Reynolds number of Re = 0.05, dimensionless shear rate of = 0.04, and power-law index of = 0.2 to 1.8. There are several interesting observations from Figure 5. First, the capsule deforms to a steady shape and then the membrane rotates around the liquid inside (tank-treading motion), which is further indicated by the streamlines in Figure 6. Second, the deformation increases with the powerlaw index. When the fluid is shear-thinning (i.e., < 1.0), the deformation is smaller compared to the Newtonian fluid case ( = 1.0), while the deformation is larger compared to the Newtonian fluid case for the shear-thickening fluid; that is, > 1.0. This can be explained by the power-law rheology. When < 1.0, the effective viscosity near the capsule is smaller compared to the Newtonian fluid, while the effective viscosity near the capsule is higher than that of Newtonian fluid for > 1.0. Based on the definition of in (16), the local is larger for > 1.0 and smaller for < 1.0 compared to that of Newtonian fluid. As presented by Sui et al. [1], a larger corresponds to a larger , that is, larger deformation of the capsule. Third, the Taylor shape parameter , which is used to quantify the deformation, increases with the powerlaw index. Finally, it is noted that y is approximately linear function of , as shown in Figure 5(b).
In order to study the Reynolds number effect on the deformation of the capsule, we simulate two additional Reynolds numbers, Re = 0.1 and 0.025. Figure 7 shows the deformation of the flexible capsule in a shear flow for Reynolds number of Re = 0.1 and 0.025, dimensionless shear rate of = 0.04, and power-law index of 0.2 ≤ ≤ 1.8. It is found that the deformation ( ) for Re = 0.1 is larger compared to the cases of Re = 0.05 and 0.025. However, the difference is quite small, implying that, in the low Reynolds number regime, for example, Re ≤ 0.1 in this work, the deformation of the capsule is not significantly affected by the Reynolds numbers used, as the inertial force is ignorable here, and the shear forces and capsule elastic forces are dominant. Therefore, the dimensionless shear rate ( ) should significantly affect the deformation of the capsule, which will be further verified by the simulations shown below by varying .
Finally, we study the shear rate effect on the deformation of the capsule by using  observations are obtained. First, the capsule deformation is sensitive to the dimensionless shear rate. This can be explained by the definition of in (16): measures the ratio of shear (viscous) forces to the stretching elastic forces, which is the dominant physical process here. A change of this ratio would cause significant difference in the capsule deformation. Second, the power-law index effect is stronger for larger , as indicated by the slopes of the functions shown in Figure 8(c). This can be explained by the fact that the physical for Re = 0.05 is shown in (c) for comparison. process changes from a shear force dominant to an elasticforce dominant process when varies from 0.1 to 0.004. For low , for example, 0.004, the elastic forces are dominant, and thus the shear force change caused by the change of the non-Newtonian rheology is smaller compared to that for larger , for example, 0.1. Finally, the deformed capsule is obviously biased from elliptical cylinder for large and ; for example, = 0.1 and ≥ 1.2. This is caused by the shear-induced torque on the deformed capsule and the decrease of effective bending resistance caused by the shear-induced elongation.
To further discuss the non-Newtonian effect, = − , which measures the local shear stress, is introduced [40]. compared to that for = 0.6. This is a further explanation of larger deformation for larger .

Conclusion
A numerical approach combining the immersed boundary method and the lattice Boltzmann method has been developed for fluid-structure interactions involving non-Newtonian fluids. Without loss of generality, the power-law fluid is taken as a representation of non-Newtonian fluids to present the method. This method couples the flexible structure (e.g., capsule) dynamics and the fluid dynamics by using the immersed boundary method and calculates the incompressible viscous power-law fluid motion by solving the lattice Boltzmann equation. In order to achieve the non-Newtonian rheology, a shear rate-dependant relaxation time is employed. The non-Newtonian flow solver has been validated by conducting a power-law flow in a straight channel. The power-law index has been varied from 0.6 to 1.4. The present numerical results show a good agreement with the analytical solutions for various values of power-law index, giving us confidence in the reliability and accuracy of the present numerical solution procedure.
To study the non-Newtonian effects on the deformation of a capsule in a power-law shear flow, we have performed simulations by varying the Reynolds number from 0.025 to 0.1, dimensionless shear rate from 0.004 to 0.1, and power-law index from 0.2 to 1.8. It is found that the capsule deformation increases with the power-law index for different Reynolds numbers and nondimensional shear rates. In addition, the Reynolds number does not have significant effect on the capsule deformation in the flow regime considered. Finally, the power-law index effect is stronger for larger dimensionless shear rate compared to smaller values.

Competing Interests
The author declares that they have no competing interests.

Introduction
The inhaled particles are the culprit for aggravating fog and haze pollution, which seriously jeopardizes the human health. The human upper airway (HUA), which consists of oral cavity, nasal cavity, pharynx, larynx, and trachea, is the main passageway for inhaled particles damaging human health, resulting in respiratory diseases range from simple common cold, snoring, and obstructive sleep apnea (OSA) to more serious pneumonia and even lung cancer [1,2]. To thoroughly understand the airflow dynamics and particles deposition pattern, it is of great significance to investigate the internal airflow characteristics within the HUA.
To date, many experimental studies have been carried out to explore airflow properties in the laboratory HUA models. Girardin et al. (1983) applied the LDA measurements to determine the velocity field in a human nasal cavity model made of cadaver at five coronal sections [3]. Heenan et al.
(2003) studied the airflow field in an idealized representation of human oropharynx (HOP) via particle image velocimetry (PIV), which clearly indicated the complex nature of HOP flow, with three-dimensionality and several regions of separation and recirculation evident [4]. Based on the hot-wire anemometry and flow visualization methods, Johnstone et al. (2004) measured the axial velocity field in the central sagittal plane of an idealized representation of human extrathoracic airway (ETA) [5]. They found that the airflow patterns in the ETA were rather complex, which involved several regions of separated, secondary, and recirculating flow. Doorly et al. The velocity distributions in both sagittal and coronal planes were obtained. The results suggested that the main stream went through the backside of larynx and trachea in inspiration and the frontal side in expiration.
Although the experimental measurements are manifestly helpful in understanding the airflow dynamics in the HUA, the lack of complete fidelity of hot wire and PIV has been clearly demonstrated [5,9]. On the other hand, the CFD is now becoming a strong tool not only to analyze the airflow patterns, but also to get the variations of temperature, pressure, and wall-shear stress in the HUA [10,11]. In recent years, the complex representative and realistic HUA models have been reconstructed. Wang et al. (2009) built the HUA based on the CT images of a healthy volunteer to quantitatively study the relationship between airflow structure and function [12]. It was observed that the air mainly passed through the floor of nasal cavity in the common, middle, and inferior nasal meatus, and the higher airway resistance and wall-shear stress were found on the posterior nasal valve. Mihaescu et al. (2008) applied the LES to depict the transitional/turbulent airflow field in the human realistic pharyngeal airway with obstruction, which revealed that the LES can provide an increased level of detail and accuracy for unsteady, separated, or vortical turbulent flow situations [13]. Likewise, Lu et al. (2014) calculated the pressure and shear stress distribution in the realistic HUA with OSA before and after treatment by the LES [14]. There existed a significant pressure and shear stress dropping region near the soft palate before treatment, and the flow resistance in the upper airway was decreased and the wall-shear stress value was significantly reduced after treatment. Recently, Wang and Elghobashi (2014) developed a direct numerical simulation-lattice Boltzmann method (DNS-LBM) to investigate the airflow in two patient-specific HUAs reconstructed from CT-scan data under inspiration and expiration [15]. The numerical results showed that the DNS-LBM solver can be used to obtain accurate airflow details and it was a powerful tool to locate the obstruction in the HUA.
In the previous studies, the inspiration and expiration behaviors in both experiments and simulations were mostly separately performed, resulting in discontinuities of inspiration and expiration in respiratory cycle. The respiratory intensity was also mostly assumed to be constant, which could also import discrepancy with the real respiratory process. Although the airflow patterns in the HUA considered of time-dependent respiratory intensity have been reported [12,16], the airflow dynamics in the HUA is still far from being completely understood. Therefore, we aim at numerically studying airflow dynamics in the realistic obstructed HUA under continuous inspiration and expiration with varied respiratory intensities. The present study would provide a theoretical guidance for the treatment of respiratory disease with obstruction in the HUA.

Problem Description.
Based on the medical CT-scan images of a Chinese male patient (38 years, BMI 25.7), the 3D anatomically accurate HUA model with obstruction resulting from pharyngeal collapse was reconstructed using the image processing software Mimics, as illustrated in Figure 1. The images were obtained in the axial plane with a resolution of 0.7 × 0.7 mm 2 , and the slice thickness is 0.625 mm. To implement the inlet and outlet boundary condition, two extensions were added on the nostrils, and another extension was added on the downstream region of larynx. During inspiration, the extensions on two nostrils were set as the inlet velocity boundaries, and the extension on the downstream region of larynx was set to be the outlet pressure boundary. During expiration, the boundary conditions were contrary to the inspiration process.
To design continuous inspiration and expiration, the numerical results of airflow field at the end of inspiration were imposed to be the initial state of expiration, and the simulation outcomes of airflow field at the end of expiration were designed to be the next initial state of inspiration. In our simulations, the variation of the inlet velocity with respect to the time was set to be periodic sine function; namely, where is the respiratory velocity, is the respiratory intensity, is the time of a complete respiratory cycle, and is the inlet cross-sectional area, which represents the area of nostril in inspiration and the area of trachea in expiration.

Numerical Method.
The software of CEM-CFD was applied to mesh the obstructed HUA model. Due to the complexity of respiratory structure, the meshes were generated by adopting the tetrahedral elements. To accurately calculate the airflow properties, a refined mesh was employed near the pharynx and larynx. The calculation results of velocity were mesh-convergent to within a prescribed tolerance (∼0.2%).
After meshing, the CFD software package Fluent 15.0 was applied to solve the flow governing equations for modeling airflow in the suffering HUA. The time step of calculation was 0.001 s, and the number of time steps was 2500 for individual inspiratory and expiratory airflow. The inlet velocity was designed to be periodic sine function, which has been defined in (1), by the user-defined function (UDF). Since the LES is a validated method for capturing transitional/turbulent unsteady, separated, or vortical flows accuracy [17], consequently it was applied to reveal such relevant flow features in the flow separation region located near the minimum crosssectional area of the airway and the downstream region. In the LES approach, the filtering operation for a variable ( ) can be expressed as where is the volume of a computational cell and the filter function ( , ) is defined as The filtering process effectively filters out eddies whose scales are smaller than the filter width or grid spacing, and the filtered Navier-Stokes equations for mass and momentum conservation can be written as where is the filtered velocity, is the filtered pressure, is the time, is the fluid density, and eff is the effective viscosity. The wall-adaption local eddy-viscosity (WALE) model was selected as the subgrid scale (SGS) model for modeling the effects of unsolvable scale of turbulence on solvability turbulence. The coupling between the pressure and velocity used the scheme of PISO, owing to its distinctive advantages in solving the transient problems. The userdefined inlet velocity was specified normal to the boundary planes of nostrils, and the static pressure was set to be zero at the outlet. In addition, the no-slip boundary condition was implemented on all the solid walls.

Experimental Validation.
Firstly, the comparison between the measured and calculated results in the center-line of cross-section was performed to validate the reliability of present CFD method. In the experiments, the velocity profiles inside HUA model were measured by the LDA. The measuring cross-sections (BC1-BC5) were indicated in Figure 2(a). Owing to the narrowing of obstructed HUA, only one measuring center-line in the cross-section BC4, which was shown in dash line, was chosen for comparison with the simulation results. The measurement was started from the posterior side to the anterior side of HUA with horizontal movement. In the simulations, the airflow field was solved by the software of CEM-CFD, which simultaneously adopted LES turbulent model to reveal the transitional/turbulent unsteady, separated, or vortical flows. The respiratory intensity in experiments and simulations was set to be constant with = 16.8 L/min, which corresponded to a constant inlet velocity of = 0.8 m/s that was calculated according to the respiratory intensity and the nostril cross-sectional area. Figure 2(b) shows the measured velocity profiles in the obstructed HUA model. The axial velocity distributions in BC1-BC5 experience a similar trend, which both is higher near the posterior side due to the strong effect of "pharyngeal jet" and becomes weak and negative near the anterior wall owing to the reversed flow effect. The maximal velocity in five cross-sections BC1-BC5 can reach approximate 5.1∼ 5.6 m/s. Figure 2(c) illustrates the comparison between the measured and calculated data in the center-line of crosssection BC4, where the radial distance has been normalized. In general, the agreement of velocity profiles between them is good. Although the velocity very near to the wall cannot be captured due to the limitation of LDA technique, the velocity profiles and values are almost the same as the calculation in entire measurable region, which identified that the present CFD method can be used to study airflow dynamics in the obstructed HUA model.

Grid Independence Test.
To verify the grid independence, three different gird sizes (1862545, 1283535, and 834353 cells) were chosen to conduct the calculations under the respiration intensity of 16.8 L/min at = 1.25 s, in which the maximal inspiration occurs in the respiration process. As shown in Figure 3, nine characteristic cross-sections (O1-O9), which represent the typical locations of nasal valve, anterior turbinates, posterior turbinates, oropharynx, pharynx, tip of epiglottis, middle of epiglottis, laryngeal, and trachea, and six characteristic points (A-E) were chosen in the obstructed HUA model. Owing to the narrowing of HUA, the measuring center-line in the characteristic cross-section O7, which is shown in dash line in Figure 3, was adopted to execute the grid independence test. From Figure 4, one can find that the agreement of axial velocity for three gird sizes is excellent from the posterior wall to the center of jet region. However, the axial velocity of the grid size of 834353 cells is much smaller than that of the other two grid sizes from the center of jet region to jet margin. From the jet margin to anterior wall, the axial velocity has little difference for all three grid sizes. In general, the discrepancy of axial velocity between 1283535 and 1862545 cells is much smaller comparing to that of 834353 cells. It seems that 1283535 cells can satisfy the request of calculations and was chosen in the rest of simulations.

Airflow Dynamics in the Obstructed HUA.
In the simulations, the respiration intensities of = 16.8 L/min (slight breathing), 30 L/min (moderate breathing), and 60 L/min (severe breathing) were conducted under continuous inspiration and expiration conditions. After two respiratory cycles, the velocity and pressure fields in circular breathing were obtained.  at these characteristic points. It can be noted that the largest amplitudes of instantaneous velocity appear at point C, which is located in the characteristic cross-section O5, in both inspiration and expiration processes. It is about 6.20 m/s at C in inspiration, slightly smaller than that of 6.40 m/s in expiration. The static pressures at six characteristic points are shown in Figure 6(b). In inspiration, the static pressure decreases from A to C, and it reaches a minimum of −6.35 Pa at the narrowest cross-section O5, where the airflow attains its maximal velocity. The changing static pressure from C to E usually creates a recirculation zone (flow reversal) due to the adverse pressure gradient. The expiration process causes occurrence of opposite phenomenon compared to the inspiration process except for the sharp drops of static pressure at C. During expiration, the minimal static pressure is nearly −11.16 Pa, much lower than that in inspiration. Apparently, there are two large pressure differences appearing before and after a flow constriction at C, indicating that the obstruction arises at pharynx, which also verifies the approach to locating the airway collapse using pressure difference between the sensors [18]. from laminar to turbulent flow there. C is located at the crosssection of pharynx; the greatest velocity occurs there with less fluctuation. The reduction of cross-sectional area between the nasopharynx and oropharynx generates a jet downstream of the restriction. Since D is still inside the region of jet flow, the velocity distribution experiences less fluctuation there. E is located in the area of jet reduction and back flow expansion; the fluctuation of the velocity gradually increases because of the strong turbulent and back flow there. With the constriction of larynx, the intensity of turbulent rapidly strengthens; there also exists substantially fluctuation of the velocity at F resulting from the violent low Reynolds turbulent and reversed flow. It can be confirmed that the airflow in the pharynx, larynx, and trachea is turbulent and unsteady. During expiration, due to the jet flow and the complication of epiglottis, the monitoring points generate varying degrees of fluctuation with the exception of F, revealing that the degree of monolithic turbulence in expiration is larger than that in inspiration.
The process of human respiration depends on the expansion and contraction of pleural, and the static pressure components at characteristic points are presented in Figures 7(b) and 7(d). Since the pressure in the nasal cavity is less than that of the atmosphere, A and B exhibit positive pressures in inspiration. Due to the contraction of pharynx, the static pressure at C turns negative and varies between 0 and −10 Pa. With the widening of airway, the static pressure from D to F gradually arises and becomes positive. Correspondingly, a recirculation zone is produced with the negative pressure drop from pharynx to downstream region. In expiration, the static pressures at D-F are almost coincident. Likewise, the highest negative pressure appears at C, and the static pressure fluctuates and tends to 0 at both A and B. Meanwhile, a similar recirculation zone is created from pharynx to turbinates.
Owing to the particular contraction of pharynx, the variation of instantaneous velocity and static pressure at C are chosen and calculated under three typical respiratory intensities, namely, = 16.8 L/min (slight breathing), 30 L/ min (moderate breathing), and 60 L/min (severe breathing), which are displayed in Figure 8. During inspiration, the trends of the velocity are similar, whose profiles all shape sinusoidal curves, although their amplitudes are very different. It can be observed that the velocities arise rapidly as the increase of respiratory intensity. During expiration, the volatilities of velocity become more obvious with the raise of respiration intensities due to the jet flow effect. The amplitudes of velocity in same respiratory intensity are approximate in both inspiration and respiration. From

Conclusion
In this work, the flow characteristics in the realistic suffering HUA resulted from pharyngeal collapse were numerically investigated based on the CT-CFD approach. The 3D anatomically accurate HUA model was reconstructed from CT-scan images of a Chinese male patient. The CFD solver with the LES method was applied to simulate the laminar-transitionalturbulent airflow in the obstructed HUA model. In addition, the LDA measurements for airflow field in the suffering HUA model were executed to validate the reliability of LES approach. In the calculations, the continuous circular breathing was accomplished by altering the boundary conditions, and the airflow field and pressure field were computed on the basis of real circular inspiration and expiration processes. The instantaneous velocity and static pressure in several selected characteristic cross-sections and points are discussed at slight breathing, moderate breathing, and severe breathing. The results reveal that airflow in the HUA model is an unsteady transitional/turbulent flow with low Reynolds number. The strong flow injection phenomenon caused by narrowing of pharynx is found in both inspiration and expiration. The strong flow separation and back flow inside HUA are also observed due to the vigorous jet flow effect in the pharynx.
In addition, it is found that the respiratory intensity has great effect on the airflow dynamics. The analysis for airflow field in the suffering HUA model could come up with great assistance for the treatment of obstructive respiratory disease. What is worth mentioning is that the present CT-CFD based computational method is applicable but not limited in modeling the flows in the circulatory and respiratory systems; it can also be applied to study the flows in the circulatory and reproductive systems, that is, sperm swimming in fertilization.

Introduction
Sorting various categories of particles from the mixture to achieve pure sample is of great importance in biological and medical engineering. With the rapid development of micro total analysis systems, small sample volume, high throughput sample processing, high efficiency, and precise particle fractionation are several representative requirements to guide the design of sorting scheme [1]. And correspondingly, a host of particle sorting techniques have been developed in these years: for example, the fluorescence-activated cell sorting [2][3][4], magnetic-activated cell sorting [5][6][7], dielectrophoresis sorting [8,9], and size-dependent sorting [10][11][12]. The last one has received a remarkable attention attributing to its promising advantages of low cost, high efficiency, and being label-free. There are four typical size-dependent sorting methods that are generally reported, the deterministic lateral displacement [10,13], the pinched flow fractionation (PFF) [14][15][16], the cross-flow filtering [17], and the inertial focusing sorting [18]. PFF is relatively simple because there is no extra and specific microstructure needed in the channel, and it has been used to sort polymer beads [14], microparticles [19], and emulsion droplets [20] and for blood cells [21] in recent years. In these above researches, an asymmetric pinched flow fractionation scheme (AsPFF) proposed experimentally first by Takagi et al. [19] is reported to perform a continuous separation and collection for 1.5∼5 m particles; it bettered the traditional PFF remarkably, while there are still some aspects that could be improved, for example, to perform a hydrodynamic analysis and further develop an active and controllable cell or particle sorter. In the present study, a numerical AsPFF cell sorter model is established with an immersed boundary-lattice Boltzmann method (IB-LBM), where the channel structure, the flow, the multiple sizes of cells, and their interactions are considered. Based on the model, cells with a prescribed size can be manipulated to enter a desired D-branch simply by regulating the flux of one D-branch (or the pressure of one outlet). The numerical results demonstrate that the numerical cell sorter is effective to perform an active and controllable cell sorting, which suggests an improved scheme of AsPFF and is valuable for guiding the experimental design of cell sorter on microfluidic chips.

Mathematical Models.
In the numerical model, the fluid motion is solved by LBM with D2Q9 lattice model. The discrete lattice Boltzmann equation of a single relaxation time model is [26][27][28] (x + e Δ , + Δ ) − (x, ) where (x, ) is the distribution function for particles of velocity e at position x and time , Δ is the time step, eq (x, ) is the equilibrium distribution function, is the nondimensional relaxation time, and is the body force term. In the two-dimensional nine-speed (D2Q9) model [29], e are given as follows: where ℎ is the lattice spacing. In (1), eq and are calculated by [26,30] eq = [1 + e ⋅ u 2 + uu : (e e − 2 I) where are the weights defined by 0 = 4/9, = 1/9 for = 1 to 4, and = 1/36 for = 5 to 8, u is the velocity of the fluid, is the speed of sound defined by = ℎ/ √ 3Δ , and f is the body force acting on the fluid. The relaxation time related to the kinematic viscosity of the fluid is in terms of ] = ( − 0.5) 2 Δ . (4) Once the particle density distribution is known, the macroscopical quantities, including the fluid density, velocity, and pressure, are then computed from = ∑ , u = ∑ e + 0.5fΔ , Although the lattice Boltzmann method is original from a microscopic description of the fluid behavior, the macroscopic continuity (6) and momentum equations (7) can be recovered from it through the Chapman-Enskog multiscale analysis [31]. Then the LBM maybe can be viewed as a way of solving the macroscopic Navier-Stokes equations: For the IB-LBM frame, the fluid motion is first solved by LBM; then the position of immersed boundary can be updated within one-time step of Δ through [32] U ( , ) = ∫ Ω u (x, ) (x − X ( , )) x, where X( , ) is the position of the cell membrane at time . U( , ) is the membrane velocity and u(x, ) is the fluid velocity. x is the lattice side length; Ω is the nearby area of the membrane defined by the Delta function (x − X) [33][34][35]: where In (9) and (10), denotes the total dimension of the model. The fluid-structure-interaction is enforced by the following equation [27,32,33,36]: where F( , ) is Lagrangian force acting on the ambient fluid by the cell membrane. In the present study, the cell model is proposed as where F is the tensile force, F is the bending force, F is the normal force on the membrane which controls the cell incompressibility, and F is the membrane-wall extrusion acting on the cell. The four force components are [33,[37][38][39] where , , , and are the constant coefficients for the corresponding force components. In (15), is the evolving cell area, 0 is the reference cell area, and n is unit normal vector pointing to fluid. In (16), X is the position of the vessel wall, and is the cut-off distance of the effective scope in the membrane-wall interaction.

Physical Model and Simulation
Setup. The geometry model of for cell sorting is illustrated in Figure 1 [19], where is the flux of a D-branch, Δ is the pressure difference between the buffer center and the outlet, and is the flow resistance produced by the microchannel. In order to allocate the flow averagely for all the D-branches under the same pressure boundary conditions, s in all Dbranches should be equal. A way to make be equal is described as two steps. First, set the pressure of all outlet to be the same. Second, change the length of the folded part of D-branches 2 and 3 until the stable flows of all outlets are equal. When sorting different size of cells, set the pressure of outlets 1, 2, and 3 to be the same, while the pressure of outlet 4 is regulatable, and the flows of D-branches can be reallocated by altering the outlet pressure. To quantify the the capacity of the reallocation of flow by regulating the flow of outlet 4, we define = out4 /(∑ 4 =1 out ), where bigger means bigger flow through outlet 4 and smaller flow through 1, 2, and 3. In addition, since the flow resistance in each Dbranch is the same, the flow is in proportion to Δ ; that is, regulation of flow can be simply realized by regulating the pressure difference; this means that also can be defined as Δ 4 /(∑ 4 =1 Δ ).

Results and Discussion
3.1. Validation. The method and model are validated carefully here by performing a simulation of flow past a stationary circular cylinder. This simulation is carried out by employing IB-LBM model. The computational domain is shown in Figure 2. The length and width of the computational domain are 1000 and 800, respectively. The center point of cylinder is located at = 301 and = 401 and the diameter of cylinder = 40. The cylinder is discretized into a series of points, and the spacing between two adjacent points is 0.6. The cylinder is handled by utilizing immersed boundary method (IB), and the feedback-force principle is adopted to compute the force density on the cylinder, which is described as [22,40] F where F(x , ) denotes the interaction force between the fluid and the immersed boundary (cylinder), 1 and 2 are large negative free constants, u(x , ) is the fluid velocity obtained by interpolation at the IB, and U(x , ) is the velocity of the cylinder expressed by U(x , ) = x / . Here, U(x , ) equals 0 because cylinder is stationary. In this case, the ratio of length of the recirculation zone and cylinder diameter , the drag force coefficient (18), the lift force coefficient (19), and the Strouhal number are calculated at Reynolds numbers 40 and 100: The results are shown in Table 1. As shown in Table 1, the present results show close agreements with the general results reported by other literatures. This means the IB-LBM model adopted in present paper is accurate enough.

Determination of the Inlet Flow Ratio .
In order to actualize the pinched flow to sort cells, it is necessary to establish an appropriate pinched segment in the transitional channel, which is able to lead all cells to move along with the lower sidewall of the transitional channel. There are three aspects for establishing the pinched segment. First, the width      Figure 3.
As shown in Figure 3, the cell center positions when leaving the pinched segment drop with the increase of , and finally they reach a relatively steady state when > 6. Although = 8 and = 10 seem to be much better, this means much higher shear stress, which may do damage to the cells. Therefore, = 6 is the choice for the present study. A D-branch choice for a rigid circular particle can be predicted by the following experimental equations [19]:

Effect of and Cell Size on D-Branch
where 0 is the width of pinched segment as marked in Figure 1, is the outflow ratio at outlet 4, is the total number of outlets, and is the particle diameter. According to the above two equations, the particle will enter the th ( = 1, 2, 3, 4) D-branch if ranges in the scope which can be described with (20) or (21), where (20) is for = 1, 2, or 3, and (21) is only for = 4.
The predicted and numerical results of the choice of D-branch which is related to the cell diameter and are exhibited in Figure 4. In these results, 11 numerical results out of 68 are found not to be consistent to the predicted results, which generally occur at the transition where the cell has approximate probability to enter two neighbouring Dbranches. A most possible reason to result in the 11 differences is the predicted results are for rigid particles while cells are flexible.
According to the results, by regulating , the 8 m and 16 m cells can be sent into any one of all four Dbranches, and some snapshots of the D-branch choice of 16 m cell are displayed in Figure 5. By contrast, the 20 m and 24 m cells can select one of three D-branches labeled 2, 3, and 4, and the 20 m cell snapshots are shown in Figure 6. The results indicate that, by simply regulating the flux of one D-branch, cells with the diameters ranging

Summary and Conclusion
A size-dependent cell sorting model with an asymmetric pinched flow is investigated numerically by immersed boundary-lattice Boltzmann method. In the model, three aspects are summarized as the following. First, the geometry of the channels is designed specially according to the effective cell sorting, where the size of the transitional channel for controlling the pinched segment is discussed in detail. Second, the parameters and are defined, respectively, for the flux ratio of the two inlets and the flux proportion of outlet 4 in all outlets. = 6 is considered as a proper value to prepare for the cell sorting, based on which the regulation of can manipulate cells with different diameters to enter different D-branches. Finally, four sizes of cells are taken into account to exhibit the capacity of cell sorting, and the relations of the regulation flux, the cell size, and the choice of D-branch are analyzed systematically. The simulation results indicate that cells with different diameters can be successfully sorted into different D-branches, this evinces that the model we established is effective, which can provide a directive reference for the design of microfluidic chip for sorting multiple sizes of cells or particles.

Introduction
The ion permeability of red blood cell membrane is associated with a variety of life activities. The hydrophilic nanopores in a cell membrane can serve as a pathway for inserting biological molecule into the cell [1,2]. The detailed understanding of ion distribution and transport in nanopores is essential to improve the understanding and controlling of ion permeability and cell transfection [3][4][5]. In another way, as microfabrication techniques continue to be developed, more and more micro/nanofluidic devices have been devised. Nanofluidic devices, such as organic or inorganic pores and channels, have primary dimensions comparable to the Debye length, and so they have been wildly used for the successful separation of DNA or biomolecular sensing down to the single-molecule level [6][7][8][9][10]. In these devices a modulation in a baseline ion current can be observed when DNA or a biomolecule is translocated through the nanopore. Analysis of the ion current modulation can be used to gather information about the specific DNA or biomolecule of interest [11][12][13][14]. As could be expected, a detailed understanding of the device ion distribution is essential to the analysis of the ionic current signals collected during nanopore-based biosensing [12,13,[15][16][17][18]. However, an in-depth understanding of the fundamental physics of ion and biomolecule behavior in the highly confined nanoenvironment of a nanosensor is far from complete. For example, a clear picture describing the complex interactions between the mobile ions in the solution, the surface charges, and the charges on the biomolecules themselves has yet to be put forth. Previously, only simple models have been proposed to explain the current modulation. The lack of accurate models to describe the transport laws of ions and biomolecules confined in nanofluidic channels not only has restricted the precision of the nanofluidic devices, but has also blocked them from more extensive application. Molecular dynamics (MD) simulations are a useful tool to study nanoscale fluid flow. By modeling and solving complex motion equations, the space, position, and velocity of each particle in the system can be defined. As a further step, the macroscopic quantities such as ion radial distribution, degree of screening, and potential profile can be analyzed quantitatively, providing a level of detail very difficult to arrive at experimentally. In this work, an MD model of a cylindrical nanopore 3 nm in radius was built and used to study the influence of solution concentration on the ion radial distribution, screening effects, and the potential profile of sodium chlorine solution confined in the pore. Simulation results indicated that as the solution concentration increased, the density peaks of both coion and counterion concentrations increased at different speeds as screening effects appeared. Due to the negative surface charges, the potential of the solution is negative near the charged nanopore wall but quickly becomes positive as the distance from the wall increases. Results from this simulation can be used to modify the current hydrodynamic model based on continuum theories and build an accurate mathematical model that can be used to describe the transport rules of ions and biomolecules confined in nanofluidic pores.

Details of the Molecular Dynamics Model
A molecular dynamics model of bulk-nanopore-bulk, which is similar to a nanopore in a cell membrane, as shown in Figure 1, was modeled for different concentrations of solution using a modified TINKER 4.2 [19] MD package. The nanopore was filled with NaCl solution, with the counterions and coions randomly distributed in the solution. The initial setting of the number of coions and counterions gave the model electrical neutrality [20]. The wall of the nanopore, however, was distributed with elementary charges along the direction which remained frozen to their original locations during the simulation [21]. The model and simulations included solution concentrations of 0.6 M, 1.3 M, and 2 M, with 1045 total water molecules used in the model. The Lennard-Jones (LJ) potential was used to approximate the interaction between a pair of atoms [20,22]. The electrostatic interactions among surface charges, ions, and water molecules were modeled using the Ewald summation algorithm [22]. The water molecules themselves were modeled using SPC/E (extended simple point charge) [23]. Table 1 gives a complete list of the parameters used for the Lennard-Jones interaction in the calculation [19,20]. The first 4 ns of the simulation were used to equilibrate the system, while the following 4 ns were used to obtain statistical data across the various solution concentrations. Qualitatively, the wall of the nanopore is similar to a hydrophilic surface. The Steele potential is used to describe the interaction between the solid surface and the fluid molecules. The values recorded for interactions between the solid surface and the molecules within the nanopore agree well with previous results reported in the literature [20]. The Steele potential of the simulation is as follows: where Δ = 2.709Å, =42.76 nm −3 , and and are obtained from bulk silica parameters: = 3.0Å and / = 230 K. The fluid molecular parameters were seen to follow the Lorentz-Berthelot rules. is the radius of the nanopore, while is the distance of the fluid molecules from the center of the nanopore. The velocity Verlet algorithm [24] was used to integrate Newton's equations of motion over 2 fs time steps. A Berendsen thermostat [25] was used to maintain the system temperature at 298.0 K with a time constant of 0.1 ps. Other simulation details can be found in a previously published paper [19]. Figure 2 shows the concentration profiles of counterions along the radial direction for different solution concentrations. Due to the negative surface charges, Na + ions are attracted to the wall surface while also being repelled by the Steele potential at the wall. The highest density peaks accumulated 1.5Å away from the charged wall. As the solution concentration increased, the density peak increased in intensity while staying at the same position. This is because the negative surface charges and charge density, as well as the electrostatic interactions among the counterions, do not change. Figure 3 shows the distribution profile of Cl − ions. These negatively charged ions are repelled by surface charges on the nanopore wall, while also being attracted by the sodium ions which have relocated near the wall. Based on these forces, the peak density is located 3.5Å away from the charged wall in the first density valley of Na + ions. As the density of Na + ions increases, more Cl − ions are attracted and accumulated, leading to increase of both density peaks, though at different speeds. Na + ion density was seen to increase more slowly than that of Cl − ions. This behavior is due to the fact that near the position of the highest density the highly packed Na + ions strongly repulse each other and therefore limit continued accumulation. In the case of Cl − ions, the peak density increases relatively faster because it is influenced not only by the solution concentration but also by the density peak of the Na + ions.

Results of MD Simulation and Discussion
Screening of the surface charge due to the ions in solution was computed using solution concentrations of 0.6 M, 1.3 M, and 2 M. The screening factor is defined as follows [14]: where is Faraday's constant, ( ) denotes the concentration of ion at the position , and is the surface charge density.
( ) > 1 corresponds to overscreening of the surface charge. Figure 4 shows the screening factor given by different solution concentrations along the nanopore radius. As the solution concentration increases, the distance between the charged wall and the point where the surface charge is overscreened decreases, and the effect of thermal motion in the electric double layer becomes comparatively weaker as more ions are shielded from the surface charges. It is obvious that the shielding effect is influenced by solution concentration. The potential within the nanopore is related to both the surface charge density and the solution concentration. Figure 5 shows the influence of solution concentration on potential. The potential is computed using Coulomb's Law and the principle of point charge superposition. The equation can be written as in which is the conversion constant, is the total number of ions in solution and the elementary charge on the surface, is the electric quantity held by ion or its elementary charge, 0 is the reference point (taken as the pore center point here) at which the voltage is taken to zero, and is the position of ion . Because of negative surface charges the calculated potential is negative near the charged wall but quickly becomes positive as the distance away from the wall increases and Na + ions have accumulated. As higher solution concentrations were used, this phenomenon becomes more evident.

Conclusions
MD simulations were used to investigate the radial ion distributions, screening factors, and potential distribution of solutions with various concentrations confined in a cylindrical nanopore. Simulation results indicate that as the solution concentration increases the density peaks of both the counterion and coion also increase, though not linearly. Because the density of surface charges remains constant, the Na + ion density peak is only influenced by the solution concentration.
Meanwhile the Cl − ion density peak is influenced not only by the solution concentration but also by the density of Na + ions. As a result, the density peak of Na + ions increases more slowly than that of Cl − ions. As the solution concentration increases, more sodium ions accumulate near the charged wall, making screening effects more evident. The calculated potential is related to both the charge density of the nanopore and net space charges in the diffuse part of the double layer. Due to negative surface charges, the potential is negative near the charged wall and becomes more positive away from the charged wall. These results can serve as an important reference value for the theoretical study of ion permeability and cell transfection as well as the design of nanofluidic devices used for the separation of DNA and biomolecular sensing.

Introduction
A cerebral aneurysm is a disease of the vascular wall causing pathological dilatations that frequently occur at the arterial bifurcations in the circle of Willis. The most serious consequence of cerebral aneurysm is the subarachnoid hemorrhage caused by the aneurysm rupture [1][2][3][4]. Since patient with unruptured aneurysm does not show any symptoms in daily life, the detection of the intracranial aneurysm and its prevention and treatment are very difficult. The recent enhancement in the technology of the imaging diagnosis systems, such as computed tomography angiography (CTA) and magnetic resonance angiography (MRA), has facilitated the detection of the unruptured aneurysm, therefore allowing for the early treatment to prevent subarachnoid hemorrhage [5,6]. There is a correlation between the size of an aneurysm and the likelihood of its rupture; for example, the larger the aneurysm is, the more likely it is that the rupture will occur [7,8]. However, because the size and the shape of an aneurysm vary with patients, evaluating the risk of rupture based on these factors is limited [9,10]. Since the prognosis of subarachnoid hemorrhage is still difficult to predict, there has been an increase in prophylactic surgery. Current treatment options include surgical clipping of the aneurysm neck and endovascular coiling, but these procedures have the risks and complications such as stroke and aneurysm rupture. So far, aneurysm size and location and multiplicity and patient's age, family history, and medical history have been considered the predominant factor to determine whether the surgery should be conducted or not. If the diameter of an aneurysm is greater than 5 mm, a surgery is generally recommended [11,12]. According to the study of Weir et al. [10] the average size of aneurysm is 8.7 mm, in which a rupture could occur.
Although the mechanisms of aneurysm initiation, growth, and rupture have not been completely understood yet, the interaction of hemodynamic forces on the vascular wall is generally accepted as the risk factor. The process of cerebral aneurysm formation, progression, and rupture is believed to be related to the hemodynamics. Endothelial cells can have the ability to detect the fluid wall shear stress (WSS) and consequently adapt their spatial organization in the hemodynamic environment. Previous studies [8,12] have shown that a uniform shear stress field causes the endothelial cells to stretch and align in the flow direction. On the other hand, the low and oscillatory shear stress fields cause the irregular shape and the loss of the cells' orientation. Thus, the low and oscillatory patterns of the WSS promote the intimal wall thickening as well as the atherosclerosis [1,4,11]. Many in vitro and numerical observations [13][14][15] have reported the effects of WSS on the development, growth, and rupture of intracranial aneurysm.
According to the high flow theory [4], the elevated WSS can initiate the wall remodeling and weaken the wall due to an abnormal shear stress field. On the contrary, the low flow theory explains that the low blood flow in the aneurysm can result in the blood stagnation in the dome. The dysfunction of the aggregation of RBC (red blood cell), the accumulation and the adhesion of platelets, and leukocytes along the intimal surface are induced by the blood stagnation. This process may cause inflammation and leads the aneurysm wall to weaken and ultimately encounter an aneurysm rupture.
The cerebral aneurysm typically occurs at the outer curvatures, bifurcations, and branching points of the cerebral arteries. Studies indicate that the blood flows in these regions are characterized by complex flow patterns, such as the strong secondary flows and flow stagnation [1,2,16]. In particular, the distributions of WSS acting on the arterial wall periodically increase and decrease in magnitude. The larger the aspect ratio of aneurysm, the smaller the WSS in the dome. The low WSS along endothelial surface coagulates the blood cells, including the red blood cells and platelets, and leads to endothelial cell dysfunction and loss. It is predicted that the low WSS can lead to an aneurysm rupture by the endothelial cell damage, thrombus, and inflammation on the aneurysm wall.
Since the hemodynamic forces act on aneurysm, the hemodynamic factors, such as WSS, pressure on the aneurysm, and the secondary flow patterns, are significant parameters of aneurysm rupture. These factors are strongly dependent not only on the geometry of the arteries but also on the shape of aneurysm. The aim of this study is to compare the hemodynamic phenomena in the cerebral arteries before and after surgery of the aneurysm under realistic conditions. To investigate the hemodynamic characteristics in brain aneurysm, the patient-specific pre-and postsurgery cerebral arterial geometries were reconstructed from medical data and blood flow has been simulated using computational fluid dynamics (CFD).

Vascular Models Reconstruction.
All medical data in the study were acquired via DSA (Digital Subtraction Angiography, Axiom Artis, Siemens Medical System, Germany) from two patients in Chung-Ang University Hospital. Starting from the medical images, the suitable geometry models were generated in order to conduct the desired blood flow simulations. Images produced from DSA are volume data consisting of 101 slices of 2,172 × 1,860 pixels' images. MIMICS software imported these DICOM (Digital Imaging and Communications in Medicine) image files to reconstruct cerebral aneurysm models. The slice thickness of the scans is 9.11 m. After loading the image files into MIMICS, the user can see images in the axial ( plane), coronal ( ), and sagittal ( ) directions. The most important first step to convert the image files into 3D model is a segmentation process. During segmentation the user creates 3D model based on the gray values within these DICOM files using thresholding. The user can apply the "region growing" tool to remove any floating pixels in the images and to separate the bone from the blood vessel. Finally the geometry models from the medical images coming from DSA shown in Figures 1 and 2 are obtained using the image segmentation and three-dimensional model creation, as mentioned. Both cases of models 1 and 2 were ruptured aneurysms, so these aneurysms were treated by endovascular coiling (refer to red circles in Figures 1 and 2)also calling it endovascular surgery.
In model 1, the blood flowing from the M1 segment of left middle cerebral artery (MCA) divides at the left MCA bifurcation and flows into the superior division of left MCA and inferior division of left MCA. In contrast, model 2 shows the blood flowing straight along the left distal vertebral artery. The possible reasons for aneurysms rupture are aspect ratio, size, and location as well as the shape of the daughter sac. In model 2 there exists the daughter sac as indicated by the red star in Figure 2(a).

Computational Models and
Simulations. The geometries generated from MIMICS are imported in ANSYS workbench to carry out the preprocessing, processing, and postprocessing. Although the blood has a Non-Newtonian behavior, it is considered a Newtonian fluid because there are no significant differences for WSS in the large arteries between Newtonian and Non-Newtonian cases [3,15,17]. Additionally, the Non-Newtonian effect can be negligible if the Reynolds number ranges from 100 to 850 [18]. Rigid cerebral arterial walls assumption is used to simplify the computation. Thus, the mass and momentum conservation equations for an incompressible Newtonian fluid neglecting the influence of body force, such as gravity, can be written as where ⃗ , , and denote blood velocity vector, blood density (1,050 kg/m 3 ), and blood dynamic viscosity (0.0035 Pa-s), respectively. The boundary conditions were as follows: since patientspecific blood flow information was not available, a physiological velocity profile was as shown in Figure 3 at the inlet [1]. Every time step the inlet velocity can be calculated using zeroth-order interpolation method (refer to ANSYS FLUENT user's guide). The physiological outlet pressure at the iliac and the renal arteries will be different. However, we did not have any pressure information at the outlet in both the iliac and renal arteries. Thus we assumed that zero pressure at the outlet boundary was employed, and no-slip condition was applied at arterial wall surfaces. The mean inlet velocity varies between 0.35 and 0.61 m/s within 1 cardiac cycle. The Reynolds numbers at the maximal, minimal, and mean rates of blood flow were approximately 220, 124, and 154, respectively. The peak systolic flow occurred at = 0.08 s as shown in Figure 3. The Womersley number, , is defined by = ( /2)√ / . represents the angular frequency of the pulsatile flow and is the diameter of inlet. In this simulation, the Womersley number was approximately 0.823. In Figure 3, five different times ( 1 , 2 , 3 , 4 , 5 ) were selected to analyze the hemodynamic behaviors. The inflow at the inlet of middle cerebral artery accelerates early in the inlet velocity temporal waveform, 1 = 0.05 s, and reaches the maximum of 0.61 m/s at 2 = 0.08 s. After this point, the velocity magnitude begins to decrease and drops to 0.42 m/s at 3 = 0.18 s. Then, the blood flow slightly oscillates in the range between 0.38 ( 4 = 0.52 s) and 0.35 m/s ( 5 = 0.88 s).
The governing equations were solved using the commercial software package FLUENT, ANSYS 15 (ANSYS Inc., Canonsburg, PA). The QUICK scheme was used to discretize the velocity and pressure variables using the PISO (Pressure Implicit with Splitting of Operators) algorithm. The implicit time-marching first-order scheme with the time step Δ = 0.01 was used for the calculations, and the maximum iterations per time step were set to 2,000.
To establish mesh independence of our numerical results, we have tested a wide range of mesh densities. For model 1 a typical mesh used in calculating the results consists of 1.56 million element cells, while for model 2 the adapted number of mesh elements is 1.87 million element cells. At this mesh density, the deviation in WSS values relative to a finer mesh is smaller than 1%. The usual mesh size, dx, is about the range of 1.74 * 10 −4 ∼8.77 * 10 −4 mm.
To obtain the stable solutions, the computation performed four cardiac cycles and the result at the fourth cardiac  cycle was used for the analysis in order to decrease numerical errors compared to the result at the third cardiac cycle less than 1%. The convergence tolerance for the continuity and velocity residuals was set at 10 −5 . The PCs used for simulations were based on the Intel5 Core6 i7-3770K processor with 3.5 GHz clock speed and 32 GB RAM memory running on the 64-bit Windows 7. The simulation time for the computational domains used in the study was approximately 96 CPU hours.     Figure 4.

Flow Structures Comparison in the
Upon comparing the flow distributions in pre-and postsurgery model 2 in Figure 5, it can be seen that the blood smoothly flows through the left distal vertebral artery except in the aneurysm region. In presurgery model 2, the result shows that the strong helical flow is generated inside the aneurysm. For postsurgery model 2, the incoming flow has the centrifugal force due to the curvature in the inlet region and blood leads to the impact on the surgical area. Figures 6 and 7 show the results of the pressure distributions on the arterial wall in pre-and postsurgery models 1 and 2. For presurgery model 1, the highest pressure acts at the corner of the left neck of aneurysm in the left MCA bifurcation. However, the highest pressure occurred at the right side of the neck of the removed aneurysm in the inferior division of left MCA for postsurgery model 1. For model 2, the peak pressure can be seen to appear in the left region of the aneurysm neck. It is believed that the inertia increases in the region of aneurysm, while the maximum pressure acts with increasing momentum due to the strong curvature of the blood vessel.
Based on the results in models 1 and 2, we could see that the inflow jet streams in the presurgery models flew into the aneurysm sac much more strongly than those in the postsurgery models. The strong, concentrated jet streams in the presurgery models had a forceful impact on the small surfaces in the aneurysm. On the other hand, the flow patterns in the postsurgery models were quite. For example, the inflow jet streams in the postsurgery models diffused more quickly. Consequently, the jet streams hit across a wider region in the aneurysm and the impact was not as severe.

Wall Shear Stresses in the Aneurysms before and after
Surgery. Since the hemodynamic forces such as WSS, pressure, and the flow patterns that act on the aneurysm are significant factors of aneurysm rupture, the WSS variations depend on the vessel geometry and aneurysm irregularity. Seo and Byun [15] mentioned that, by increasing the aspect ratio of aneurysm, the WSS around the dome of aneurysm gets smaller and damages endothelial cells, which results in the inflammation of the wall of aneurysm. Consequently, one can predict that an aneurysm rupture will occur.
If the sizes of aneurysms are greater than 10 mm and the aspect ratio is larger than 1.6 [10], aneurysms have a high probability of rupturing. The aspect ratio, defined as (Line C-D)/(Line A-B), of an aneurysm is 2.04 in model 1, while the aspect ratio of model 2 is 1.58. According to these criteria, model 1 is likely to rupture compared to model 2. Figure 8 shows the WSS distributions along the yellow lines (Line 1) of preaneurysm and postaneurysm sides at five distinct time points. For preaneurysm model 1, the magnitude of WSS at the throat of superior division of left MCA has the peak value due to the impingement of blood flow into the aneurysm. Inside of the dome, WSS will dramatically drop due to low flow velocity and flow recirculation. The average value of WSS along the dome line has 3.64 Pa for peak flow (at 2 ), while the average WSS has 2.14 Pa at 5 . Fukazawa et al. [19] reported that the average WSS at the dome was 2.27 Pa. For postaneurysm model 1, the magnitudes of WSS increase to 43.8 Pa at 2 and 19.7 Pa at 5 . Previous studies [20,21] have shown that the magnitudes of WSS in healthy middle cerebral arteries were approximately 18 to 22 Pa. In particular, Lindekleiv et al. [22] found the maximum WSS of 33.17 Pa at the female MCA bifurcation with the blood velocity 0.75 m/s. It was even hard to compare our patient-specific models with Lindekleiv et al. idealized models [22]. The maximum WSS in the MCA was 24% higher in our simulation compared with the result of Lindekleiv et al. [22].
In contrast to model 1, the very low WSS in model 2 of less than 2.5 Pa can be seen in the dome area. The WSS is low due to the slow blood velocity in the dome region, while the blood flow in model 1 is separated in both directions in places where the blood flow momentum can impinge. However, it can be seen that the WSS acting on the blood vessel after operation for model 2 has the maximum 107.36 Pa and the minimum 42.37 Pa. For the comparison of the geometry shapes between models 1 and 2, it can be seen that the WSS in normal blood vessel of model 2 after surgery is about three times larger than that of model 1.
In our study the low WSS found in the dome along line 1 in Figure 1 is consistent with CFD study by Shojima et al. [12]. Figure 9 shows the WSS distributions for models 1 and 2 before surgery at four different times. In Figure 9 the WSS in model 1 around the aneurysms surface at peak flow rate, = 0.08, was unevenly distributed. However, the WSS distributions changed considerably with complex flow patterns (see Figure 4(d)) during the diastolic phase. Further, the whole area of aneurysms was exposed to low WSS.
Since the pressure gradient (see Figures 7(a)-7(d)) has been favorable in the left distal vertebral artery, the blood flow in the region is uniform as shown in Figures 5(a)-5(d). As a result, during the whole cardiac cycle the WSS on the left distal vertebral artery is homogeneously distributed in model 2. The heterogeneous distributions of WSS are shown in the aneurysmal regions (see Figures 9(e)-9(h)). As the blood velocity decreases, it can be seen that the regions of low WSS on the daughter saccular aneurysms which is indicated by the red star in presurgery model 2 in Figure 2 widen. It is more likely to rupture the aneurysms in this region.

Conclusions
The purpose of this study was to investigate the hemodynamic characteristics in the cerebral arteries before and after    surgery of the aneurysms under realistic flow conditions. We performed computational fluid dynamic simulations for the geometries generated from CT image-based patient-specific geometries taken before and after surgery. The two patients in the study underwent cerebral angiography with 3D reconstruction using image processing software for 3D design and modeling developed by Materialise. The differences in the unsteady flow, pressure, and WSS between pre-and postaneurysms models were characterized in this study.
Even when the blood vessel geometry was quite different, this study found that the flow patterns with the complex vertical structures inside the aneurysm were similar. The results of flow patterns for models 1 and 2 after the surgery were quite different. It was also discovered that the highest pressure occurred at the distal neck of aneurysms. The average wall shear stress after surgery for model 1 was approximately three times greater than that before surgery, while it was about twenty times greater for model 2. The comparison between models 1 and 2 showed that the WSS in normal blood vessel of model 2 after surgery was about three times larger than that of model 1.
The maximum WSS appeared at the proximal superior division of left MCA in model 1 and at the distal aneurysm neck in model 2. During the diastolic phase the magnitudes and distributions of WSS in aneurysm region in models 1 and 2 were markedly reduced and widen. In particular, the area of low WSS in the daughter saccular aneurysm region in model 2 is associated with aneurysm rupture. Thus the distribution of WSS in aneurysm region provides useful prediction for the risk of aneurysm rupture.
There are limitations to this study. Although the vessel wall is elastic in the study, we do not consider the elasticity of  (Pa) Figure 9: Wall shear stress distributions for models 1 and 2 before surgery at four different times. the blood vessel. However, it is believed that the compliance of the blood vessel may have an important role in influencing the hemodynamic changes. The boundary conditions at the inlet and outlet were not patient-specific but were applied from Takao et al. [1]. The study also only looked at two patients for the simulation with the atmospheric pressure at both outlets. With consideration of the elasticity and the real physiological conditions, we expect to perform the simulation for future analysis.