|Home | About | Journals | Submit | Contact Us | Français|
Immobilization of cells inside microfluidic devices is a promising approach for enabling studies related to drug screening and cell biology. Despite extensive studies in using grooved substrates for immobilizing cells inside channels, a systematic study of the effects of various parameters that influence cell docking and retention within grooved substrates has not been performed. We demonstrate using computational simulations that the fluid dynamic environment within microgrooves significantly varies with groove width, generating micro-circulation areas in smaller microgrooves. Wall shear stress simulation predicted that shear stresses were in opposite direction in smaller grooves (25 and 50 μm wide) in comparison to those in wider grooves (75 and 100 μm wide). To validate the simulations, cells were seeded within microfluidic devices, where microgrooves of different widths were aligned perpendicularly to the direction of the flow. Experimental results showed that, as predicted, the inversion of the local direction of shear stress within the smaller grooves resulted in alignment of cells on two opposite sides of the grooves under the same flow conditions. Also, the amplitude of shear stress within microgrooved channels significantly influenced cell retainment in the channels. Therefore, our studies suggest that microscale shear stresses greatly influence cellular docking, immobilization, and retention in fluidic systems and should be considered for the design of cell-based microdevices.
Microfluidic devices are enabling platforms that can be useful for studies related to high-throughput screening, biochemical synthesis, and biology1,2,3. These devices manipulate fluid flows, minimize costly reagent consumption, and enable high-throughput experimentation in a controlled manner4,5,6. Microfluidic devices can also be used to control the cellular microenvironment by regulating the temporal and spatial presentation of soluble factors as well as shear stress on cells7, 8,9. To immobilize cells within microfluidic channels, a number of approaches such as encapsulation within photocrosslinkable polymers10, adhesion to patterned proteins11, and protein coatings12 have been widely used. In addition, microstructures that enable the capture of cells have been used to immobilize cells within fluidic channels. For example, cup-shaped microstructures that enable individual cell capture within specific regions of a channel have been fabricated13,14. An alternative approach to capture cells within fluidic systems is through the use of grooved substrates that create regions of low shear stress15,16. Within these systems, cells that are localized within the shear protected regions remain entrapped while cells outside these regions are washed away. Microstructure based cell docking has been used for a number of applications including liver bioreactors15, drug discovery chips16 as well as cell staining devices17. Intrinsic to the ability to analyze cells is an understanding of the role of shear stress on capturing cells and its role on the behavior of immobilized cells.
To study cellular behavior under simulated physiological microenvironments, it is necessary to control conditions such as mechanical stimuli18. Shear stress induced on the cell membrane by medium perfusion in a microfluidic device can play an important role in cell growth, differentiation, and metabolism19. Shear stress generated within flow-based microfluidic devices has been shown to affect cell docking and adhesion20. Computational modeling has been extensively used to estimate the shear stress distribution within the complex microstructures of tissue engineering scaffolds and bioreactors21,22,23,24 and to evaluate the fluid behavior in microfluidic devices such as micropumps and micromixers25,26. In microchannels or microvessels with square or circular cross-section, the evaluation of the shear stress acting on the main channel walls can be analytically carried out by using Poiseuille models18,27,28. However, in these devices, the shear stress pattern in the docking sites (i.e. microgrooves or microwells) can not be determined by simple analytical calculations and requires computational approaches.
Despite an increasing number of studies that utilize microgrooved substrates to immobilize cells within microchannels15,16,17,29, relatively few studies have analyzed effects of the dynamics of fluid flow within grooved substrates30. Thus an experimental and theoretical framework to analyze the effects of the unique fluid flow properties in these regions, such as the presence of micro-circulation areas is of value for understanding cell alignment and function within such cell-based microdevices.
In this work, we generated computational simulations to predict the magnitude and direction of shear stress within microgrooved channels. Furthermore, we developed approaches to validate computational simulations on cell alignment and docking in the microgrooves. Computational Fluid Dynamics (CFD) modeling revealed the effect of various geometrical configurations (i.e. microgroove widths) and system parameters (i.e. inlet flow rates) on the fluid dynamic environment (i.e. shear stress distribution and micro-circulation areas). These results were confirmed with experiments in which the seemingly perplexing cell alignment behavior in grooved substrates is explained. Therefore, our results demonstrate that micro-circulations within microgrooves can be used to align cells in microfluidic channels and provide a theoretical and experimental framework for the design of microdevices containing grooved substrates for cell capture.
Microfluidic devices were fabricated by using previously published soft lithographic methods1, 16, 31, 32. Master molds patterned with 40 μm thick resist were made by patterning a negative photoresist (SU-8 2015, Microchem, MA) on a silicon wafer. A negative replica of microchannels in poly(dimethylsiloxane) (PDMS) (Sylgard 184 Silicon elastomer, Dow Corning, MI) was fabricated by replica molding against the master mold. Briefly, PDMS molds were generated by mixing silicone elastomer and curing agent (10:1 ratio). The PDMS prepolymer was poured on the silicon master that was patterned with photoresist and cured at 70°C overnight. PDMS molds were then peeled off from the silicon wafer. Cell inlets and outlets were punched by sharp punchers for medium perfusion and cell seeding. Figure 1 shows the schematic design of a microfluidic device. The microfluidic device consisted of two PDMS layers: a top fluidic channel and a bottom microgrooved surface. The top fluidic channel was 40 μm in height, 5 mm in length, while the bottom microgrooves, which were placed perpendicularly to the flow direction, were 40 μm in depth, 4 mm in length, and either 25, 50, 75, and 100 μm in width. An array of 67 microgrooves was localized within a microfluidic channel for each width, patterned at equal spacing as their widths. The two PDMS layers were aligned and bonded after each surface was treated with oxygen plasma (10 min at 30W, Harrick Scientific, NY).
Cardiac muscle cell line (HL-1) previously derived from the AT-1 mouse atrial cardiomyocyte tumor lineage33 was used to study cell alignment under medium flow. Cells were cultured with medium which contained 87% Claycomb Medium, 10% Fetal Bovine Serum (FBS), 1% Penicilin/Streptomycin, 1% Norepinephrine, and 1% L-Glutamine. The cells were cultured at 37°C culture in a humidified 5% CO2/95% air incubator. NIH-3T3 mouse fibroblasts were also used to study cell alignment and cell retention. They were cultured in Dulbecco's Modified Eagle's Medium (DMEM) containing 10% FBS. Tissue culture medium and serum were purchased from Gibco Invitrogen Corporation, CA. All chemicals were purchased from Sigma, Mo.
To seed the cells into the microfluidic device, the cells were trypsinized and dissociated with culture medium. The microchannel was filled with medium and the cells were seeded through a cell inlet port at the cell density of 4×106 cells mL−1 that allowed uniform cell distribution. Subsequently, medium was infused using a syringe pump (PHD 2000 Infusion, Harvard Inc, USA) at an average flow rate of 5 μL min−1, corresponding to a calculated inlet velocity of 5.2×10−4 m sec−1. Cell adhesion to the bottom of the grooves was minimized by starting the experiments within 10 minutes after the seeding. In addition, different inlet velocities (10.4×10−4, 20.8×10−4, 52.1×10−4 and 78.1×10−4 m sec−1) were used to investigate the influence of shear stress on cell retention within microgrooved channels. For the cell retention studies, cells were first seeded within the microgrooves, and then exposed to 5 μL min−1 flow to obtain stable alignments at the base of the grooves. The inlet velocity was then progressively increased to the fixed velocities presented above and the washing at each step was performed until no cells were washed out (i.e. steady state conditions). The medium was infused from smaller grooves (25 μm width) towards wider grooves (100 μm width). The opposite flow direction was also used to test the sensitivity of the results with respect to the flow direction.
Phase contrast images were taken by using an inverted microscope (Nikon TE 2000-U, Nikon Inc., USA). For each condition and groove width at least 30 grooves were analyzed. For analyzing the localization of cells to either corner of the channel, the number of cells distributed on the upstream corner of the microgrooves was compared to the number of cells distributed on the downstream corner. These data were compared with the total cell number in the microgroove to calculate the percentage of the cells at each side of the grooves. The number of the cells for each groove varied from at least 30 cells in 25 μm to at least 140 in 100 μm wide grooves. The total number of cells for each groove width was normalized to be 100 %. For the cell retention experiments, the total cell number for each microgroove width was counted for each inlet velocity, after reaching steady state conditions. The experiments were performed in triplicate for each flow condition. At least 30 grooves were analyzed for each groove width (25, 50, 75 and 100 μm), both after cell positioning (inlet velocity of 5.2×10−4 m sec−1) and in the cell retention experiments (after the application of inlet velocities of 10.4×10−4, 20.8×10−4, 52.1×10−4 and 78.1×10−4 m sec−1). Statistical analyses were carried out using the t-test with p<0.05 considered significant.
Computational fluid dynamics was used to predict the wall shear stress as a function of the channel width and the flow rate using FEMLAB software (COMSOL 3.2, Comsol AB). Since the channel length (5 mm) in the device was significantly greater than the channel height (40 μm), the system could be successfully modeled using a 2D simulation. We performed the modeling of a microfluidic device with microgrooves of different widths (25, 50, 75, and 100 μm). The culture medium was modeled as an incompressible, homogeneous, Newtonian fluid with the same properties of water, with a density (ρ) of 1,000 kg m−3 and a viscosity (μ) of 0.001 Pa s, as in previous studies34. The steady state Navier-Stokes equations for incompressible fluids were solved using COMSOL:
where v and p are the velocity vector and pressure, respectively. No-slip boundary conditions were applied for the channel and groove walls and a velocity profile was applied to the inlet boundary. Different velocity inlets (i.e. 5.2×10−4, 10.4×10−4, 20.8×10−4, 52.1×10−4, and 78.1×10−4 m sec−1) were used to simulate the experimental conditions. A zero pressure condition was applied to the outlet. The fluid domain was meshed using De Launay algorithm in COMSOL 3.2, resulting in 38052 degrees of freedoms and 8072 triangular elements. The average triangular element was 5 μm on a side with a minimum element size of 1 μm to resolve wall shear stress distribution. The mesh was progressively refined through mesh sensitivity analyses: at each simulation the elements showing high velocity gradients were refined, until reaching convergence of sensitive measures of the predicted quantities (error below 5% on shear stress values at the walls). As a result, triangular elements, 5 μm on a side, were mostly used in the center of the microchannel, where the velocity gradient was lower, and small triangular elements were used close to the walls, where the velocity gradient was higher. For validation, computationally derived shear stresses on the main channel walls were compared with the average shear stress (τfree-wall) assessed by analytical calculations based on Stokes flow and lubrication theory. In these conditions, the maximal shear stress applied to the cell (τmax-cell) was also estimated by Equation (2) simplifying the adhered cell as a hemisphere27,28.
where H and W are channel height and width, μ is medium viscosity and Q is the flow rate.
To experimentally visualize cell docking and immobilization within grooved microchannels, we developed a PDMS-based microfluidic device consisting of a fluidic channel and grooved substrates (Figure 1). PDMS molds with the impressions of the fluidic channel and grooved channels that contained different widths (i.e. 25, 50, 75, and 100 μm) were aligned and bonded by oxygen plasma treatment. Cells were seeded in the channels through medium perfusion and immobilized within the low shear stress regions in the fluidic channels by being captured in the microgrooves. In our system, the fluidic channels could be easily fabricated and visualized to enable rapid analysis of multiple experiments. Furthermore, all the grooved features were fabricated on the same array enabling us to test all the conditions within the same experimental setup.
To predict velocity profiles and shear stress patterns within grooved channels, we used computational fluid dynamics. Figures 2A and 2B represent velocity contours and streamlines of the laminar flow resulting from the 2D simulation using an inlet velocity of 5.2×10−4 m sec−1. As shown in the images, higher flow penetration was generated in large grooves (75 and 100 μm in width) compared to small grooves (25 and 50 μm in width). In addition, the streamline analysis of velocity profiles predicted the formation of micro-circulation areas within microgrooved channels (25 and 50 μm in width). Under these conditions, the direction of the local velocity near the base of these grooves was opposite to the mainstream fluid flow. Large grooves showed small micro-circulation areas present only in the corners, whilst the local velocity near the remaining surface of the groove base had the same direction as the mainstream flow. In addition to velocity profiles, we also simulated the resulting shear stresses within the microgrooves. Figure 3A represents the shear stress profile for the different widths at the base of the grooves. The inversion of the local velocity near the groove base between 50 μm to 75 μm grooves resulted in the inversion of the shear stress profile from negative to positive values. This change represents a change in the direction of the shear stress confirming the previous results from the computed velocity profiles. As expected, the magnitude of shear stress on the walls outside the grooves was significantly higher than that of shear stress within the grooved channels (Figure 3B). Thus, for an inlet velocity of 5.2×10−4 m sec−1, the wall shear stress outside the grooves was ~ 0.65 dyne cm−2, which was considerably higher than the average wall shear stress amplitude in the 25 μm wide grooves (−3.8×10−4 dyne cm−2). Furthermore, the increase of grooved widths resulted in increased average shear stresses (−1.6×10−2, 0.7×10−2, and 4.5×10−2 dyne cm−2 for 50, 75 and 100 μm grooves in widths respectively). Cells were seeded and aligned into the microgrooves using a single flow rate (5.2×10−4 m sec−1). We used this inlet velocity, since it resulted in cell alignment with low shear stresses which were not considered to induce significant changes in cell phenotype for the short durations required for the alignment. Also, Figure 3B reveals that wall shear stress values correlate linearly with the inlet velocity, consistent with the properties of low Reynolds number (0.04−0.6), which result in negligible convective terms in the Navier-Stokes equations. All the inlet velocities that were experimentally applied (5.2×10−4 m sec−1 in the cell positioning and 10.4×10−4, 20.8×10−4, 52.1×10−4 and 78.1×10−4 m sec−1 in the cell retention tests) were simulated by means of the CFD model. All the considered inlet velocities show laminar flow, as expected by the low Reynolds number applied (<1). In particular, all the studied conditions are Stokes conditions, in which the wall shear stresses change linearly with the velocity inlet, as shown in the simulations (Figure 3B).
The micro-circulation and direction of the localized shear stress generated near the groove bases can control and manipulate cell localization and alignment. Specifically, we were intersected in whether the local shear stress profiles in the channels affected cell alignment in the grooves. Figure 4A and B show the alignment of fibroblasts and HL-1 cardiomyocyte cells in 50 μm and 75 μm wide grooves. It was observed that cells were docked within the microgrooves, where the predicted shear stresses were lower than those experienced outside the grooves. Interestingly, our experimental results confirmed the predicted inversion in the direction of the shear stress at the base of the grooves that were between 50 and 75 μm wide. In the 50 μm wide grooves, cells were aligned towards the upstream corner of the groove, counter to the direction of the mainstream flow. On the other hand, within 75 μm wide grooves, cells were aligned in the downstream corner of the grooves. Figure 4C and D show time-lapse images in the 50 μm and 75 μm groove channels before reaching steady state. These time-lapse images revealed the different direction of the cell movement at the 50 μm and 75 μm groove channels. We also quantified the distribution of cells in the grooves using images from the experiment. Figure 5 shows this quantitative analysis for cell alignment on the upstream corner and downstream corner of the grooved channels. The extremely low shear stresses that were predicted within the 25 μm wide grooves did not influence cell alignment. On the other hand, the higher shear stresses within larger grooves (50, 75 and 100 μm wide) resulted in cell alignment on the upstream and downstream corners according to the localized shear stress direction and micro-circulation.
Shear stress significantly influences not only cell alignment but also cell retention within microgrooved channels. Figure 6 shows the experimental and computational analysis of cell retention as a function of inlet velocity, hence the shear stress within the channels. Cell retention in grooved channels was analyzed as a function of different inlet velocities and groove widths. This analysis revealed that increasing the inlet velocity from 10.4 to 78.1×10−4 m sec−1 significantly decreased cell retention within grooved channels irrespective of the channel width. For 78.1×10−4 m sec−1 inlet velocity, 40% of the initial number of the aligned cells was remained (Figure 6A). In addition, we also simulated the wall shear stresses for different inlet velocities for various microgroove widths (Figure 6B). As expected, the predicted shear stress increased as a function of the increasing flow rate. Thus, cell retention within grooved channels decreased with increasing shear stress and inlet velocity. Therefore, localized shear stress and micro-circulation generated within microgrooved channels controlled and influenced cell alignment and retention.
This study provides a computational and experimental platform for the study of fluid flow properties on cell alignment and retention within grooved microfluidic channels. Computational fluid dynamic modeling was used to predict cell positioning within an experimental setup of a microfluidic device. In particular, this work shows agreement between experimental cell alignment on the groove corners and magnitudes and directions of computationally predicted shear stresses within the grooves. It was demonstrated that the localized shear stress direction and micro-circulation significantly controlled cell alignment within grooved channels. In addition, shear stress profiles were influenced by the groove widths. Small grooves (25 and 50 μm in width) resulted in the formation of micro-circulation areas, which inverted the direction of the shear stress near the base of the grooves. In 50 μm wide grooves, where the shear stress was high enough in magnitude to move the cells, the cells aligned on the upstream corner of the grooves, while cells resulted in being randomly distributed in the 25 μm wide grooves where the shear stress was extremely low. In large grooves (75, 100 μm wide), there was higher penetration of the mainstream flow and small micro-circulation areas were present only near the corners. The wall shear stress at the center of the base of the grooves had the same direction of the mainstream flow, aligning the cells on the downstream corner. These results were achieved by two different cell types (i.e. cardiomyocytes and fibroblasts) with an inlet velocity of 5.2×10−4 m sec−1.
The presence of micro-circulation areas in 50 μm wide microgrooves was in agreement with previous findings from Horner et al.30, whom studied oxygen transport within a microgrooved bioreactor. In this work, the computational modeling was carried out for a single inlet velocity, comparable to the one applied in our study. However, as compared to our study, micro-circulation areas were studied only in 50 μm wide microgrooves. We found that micro-circulation areas also occur in smaller grooves (25 μm in width), but not in wider ones (75 and 100 μm in width), and build on their work to demonstrate the effects on wall shear stresses and subsequent alignment of cells in microchannels.
In our study, the computational model was validated by comparing the computed shear stresses at the channel walls, outside the grooves, with the corresponding analytical values. Indeed, when the Reynolds number is low (<< 1), the flow into channels can be approximated to a Stokes flow and the average wall shear stress (τfre-cell) is usually estimated by solving Equation (2). Considering the channel geometry and an inlet velocity of 5.2×10−4 m sec−1 (5 μL min−1 flow rate), the analytically calculated wall shear stress on the channels is 0.63 dyne cm−2, comparable to the results found by the computational model (0.65 dyne cm−2). Computationally predicted wall shear stresses within the microgrooves were found to be significantly lower than in the main channel as reported in previous studies29, in particular when employing comparable set-up conditions15. Further computational models, simulating several different groove widths within the 50−75 μm range should be carried out and eventually validated through experimental tests. This work could be addressed as further development of the present research.
One of the main limitations of our computational model is that it does not take into account the actual geometry of the cells and their distribution within the microgrooves: the wall shear stresses are considered to be an index of the shear stresses acting on the cell membrane, although the cell presence is not simulated. However, for the low magnitude of Reynolds numbers, the cell presence does not strongly change the flow. It has been already calculated by previous studies that the average shear stress on the cell membrane, when the cell adheres to a channel, is similar to the average wall shear stress without the cell27. This approximation may be also applied to the larger microgrooves (50, 75, 100 μm wide), where the flow near the base of the grooves can be approximated to a stationary flow into the channel (Poiseuille flow) and the cell size (15 μm diameter35) is relatively small compared to the groove. On the other hand, the cell size in the small grooves (25 μm wide) becomes more relevant as compared to the groove dimensions. Thus, actual shear stress experienced by the cell can not be analytically derived and was not taken into consideration by the computational model. A reference cell diameter of 15μm was considered in the simulations because it was in agreement with the average cell size experimentally observed for the fibroblasts used in the retention tests.
We used the inlet velocity of 5.2×10−4 m sec−1, since it resulted in cell alignment with low shear stresses in all the grooves (~ 0.045 dyne cm−2). These shear stress levels were not considered to induce significant changes in cell phenotype for the short durations required for the alignment (~ 10 minutes). After cell alignment, the inlet velocity should be accurately set to optimize the fluid dynamic environment and shear stress stimulation for the particular cell type. For instance, it is known that steady shear stresses in the range of 10−15 dyne cm−2 stimulate vascular endothelial cellular responses that are essential for endothelial cell function36. Different results were found for human chondrocytes, in which high shear stress (16.4 dyne cm−2) was found to down-regulate the expression for extracellular matrix production37, while low shear stress regimes (below 0.1 dyne cm−2) increase extracellular matrix synthesis in engineered cartilage constructs38.
Moreover, the inlet velocity should be optimized to avoid cells from being washed off in the device. After cell alignment, cell retention was tested in response to different flow rates. The cells positioned within the grooves were progressively washed off with increasing the shear stresses. The cell retention in the grooves was significantly influenced by the magnitude of shear stresses near the walls (where cells ware located) but not by the groove widths. To explain this finding, the average shear stress values experienced by the aligned cells near the groove corners (~ 15 μm from the corners) was calculated in 50, 75, and 100 μm wide grooves for the tested inlet velocities (Figure 6B). For the 25 μm wide grooves, the average shear stress was calculated on the entire base of the grooves, because cells were not significantly aligned on the corners. For the same inlet velocities, the amplitude of the considered shear stresses were similar for the 50, 75, and 100 μm wide grooves (0.1 ± 0.01 dyne cm−2 at 10.4×10−4 m sec−1), but lower for the 25 μm wide groove (7.6×10−4 dyne cm−2 at 10.4×10−4 m sec−1). This could be explained by the exclusion of cells from the CFD simulations. As expected, an increase in the calculated shear stresses as a result of the increasing flow rate resulted in increased shear stresses and lower cell retainment in the channels.
It is interesting to note that previous work18 that focused on the effect of the shear stress on cell retention have found higher percentage of cell adhesion (> 80%) after 10 minutes, when the cells were exposed to higher shear stresses than those applied in the present study. This difference is likely caused by the different levels of cell adhesion to the surfaces. In the present study, the flow was applied to channels after 10 minutes of cell seeding, not allowing cell attachment on the surface. On the contrary, in the previous study, cells were adhered on fibronectin coated microchannels.
A potential limitation of the cell retention test was the device design, since the microchannels had grooves with different widths in series (25~100 μm), allowing for possible cross-contamination. Thus, theoretically, a cell washed by a small groove may be trapped into the larger ones. However, additional experiments performed by inverting the direction of the medium flow were not significantly different from those presented in Fig. 6A (Data not shown). Cell positioning changed in the sense that the cells aligned in the opposite corner as compared with the experiments performed with the reference flow direction, as expected (i.e. counter flow in 50 μm wide grooves and in the same direction of the flow for 75 and 100 μm wide grooves).
This simple integrated computational and experimental approach can be a powerful tool for the design of microfluidic devices with controlled fluid dynamic environments. Similar integrated approaches are currently applied for the development of microdevices for different applications such as polymerase chain reaction (PCR) analyses39 or dynamic cultures in bioreactors40,41. This paper shows the feasibility of designing microfluidic devices which enables controlled cell alignment. For example, the predicted fluid dynamic fields and subsequent cell patterning within microgrooves may be useful for studying co-cultures. The control of co-cultures of different cell types can be performed by using a simple setup (i.e. different groove channels and flow direction of medium perfusion) as compared to previous methods42,43,44,45. Different cell types could be aligned on opposite corners of the same microgrooves by controlling flow direction and groove width. Therefore, computational and experimental approaches for controlling localized shear stress and micro-circulation could be useful tools to understand cellular docking and alignment in a microfluidic device.
We demonstrated, by using both computational and experimental approaches, that local shear stresses influence cell localization and retention within microfluidic channels containing grooved substrates. The amplitude and direction of the shear stress in microgrooves can be manipulated by geometry configurations (i.e. groove width) and system parameters (i.e. inlet velocity). Furthermore, we demonstrated the existence of localized micro-circulation regions inside microgrooves that appear in a size dependent manner and regulate cell alignment within grooved substrates. Also, cell retention within the grooved substrates was dependent on fluid flow rate indicating the importance of this factor in the design of grooved substrates in fluidic devices. Thus, our results provide a framework for understanding controlled patterning and docking of cells within microgrooved channels as a powerful tool for studying cell-based diagnostics, high-throughput screening, and tissue engineering.
This work was partly supported by the Coulter Foundation, National Institutes of Health (NIH), the Center for Integration of Medicine and Innovative Technology (CIMIT), US Army Core of Engineer and the Charles Stark Draper Laboratory. S. Shrivastava was supported by the Ignited Mind Fellowship Award. A. Manbachi was supported by the Canadian Natural Sciences and Engineering Research Council USRA. M. Cioffi and M. Moretti were supported by the Progetto Roberto Rocca Collaboration. The authors thank Professor Edward Hæggström, Dr. Hossein Hosseinkhani, Dr. Mohsen Hosseinkhani and Professor Gabriele Dubini for helpful discussion. We also thank Dr. Satoshi Jinno for help in illustration of Figure 1.