Electrochemical and Biological Simulation for Microfluidics

My research as a graduate student focuses on biosensing applications using electrochemical impedimetric methods. Unlike mechatronic systems, these applications consider a dynamic environment at such a microscale, it is quite hard to perceive what is really going on when we merely recognize the change of physical property. I usually ask myself: What is the underlying mechanism? Tucked away from the limit of horizon of the human eye, we often see nothing happening when doing these experiments.

That is why simulation is so important! In this project, I performed 3 simulation tasks in order to visualize the changes of important physical properties for my study (impedimetric microfluidic chip for biosensing) using COMSOL Multiphysics.

The tasks are:

  1. Velocity field simulation inside microfluidic channel.
  2. Real-time molecule immobilization on gold surface.
  3. Electrochemical impedance spectroscopy (EIS) simulation.

The result of task 1 supports the hypotheses of task 2 and 3.

Physical Environment Setting

The simulation environment is set as the interior of a microfluidic channel with gold microelectrodes. See the film below for visualization.

The chip is fabricated using soft lithography and photolithography. The microfluidic channel has a width of 1mm and a height of 100μm at the center. The gold microelectrodes form a pair for square pads (300×300μm) at the center of the channel (Fig. 1).

Figure 1. Microfabrication process, dimensions, and microscopic view of the microfluidic electrode chip used in this project.

Velocity Field Simulation Inside Microfluidic Channel

The objective of this task is to simulate the fluid velocity field inside the channel on a sliced plane. Due to the fact that 3D simulation is time-consuming, if a 3D environment can be reduced to a 2D environment, then a large amount of time can be saved. By performing this task, it can be seen if dimension reduction modeling of task 2 and 3 are feasible.

Considering the physical nature of the microfluidic channel, a laminar flow model is implemented along with the Navier-Stokes equation:

$$\rho(\textbf{u}\cdot\nabla)\textbf{u}=\nabla\cdot[-pI+\mu(\nabla \textbf{u}+(\nabla \textbf{u})^T)]+F$$

, which means a balance between inertia (\(\rho(\textbf{u}\cdot\nabla)\textbf{u}\)), pressure (\(-pI\)), viscous (\(\mu(\nabla \textbf{u}+(\nabla \textbf{u})^T)\)) and external (\(F\)) forces. A stationary study is implemented (\(\rho \nabla \cdot \textbf{u}=0\)), and water is set as the fluid.

Fig. 2 shows the 3D velocity field inside the channel. A bigger arrow indicates a larger velocity magnitude. The reason that the velocity is faster at the center is because of the lower channel height.

Figure 2. 3D velocity field inside microfluidic channel.

Fig. 3 shows animated 2D velocity fields of sliced planes inside the channel. Due to laminar flow, velocities near the boundary get close to zero. However, the steady-state velocity reaches a constant value away from the boundary.

Figure 3. 2D velocity fields for xy and zy sliced planes.

A top view of the velocity field and the velocity at different x positions are shown in Fig. 4. It can be concluded that at the center of the channel where the microelectrodes lie, the fluid velocity stays constant, and subsequent tasks can be carried out using 2D models.

Figure 4. Top-view and x position-dependent velocity magnitude of fluid.

Real-time Molecule Immobilization on Gold Surface.

For microfluidic electrochemical biosensors, it is quite common that immobilization of sensing elements takes place at the center of the channel on electrode surfaces (e.g. Au). In this task, molecules are simulated to flow past the channel, and be immobilized on the electrode pad at the bottom center. A slice on the zy plane is used as the modeling geometry of this task (Fig. 5).

Figure 5. Modeling geometry used in task 2. A sliced area of the microfluidic channel with pad electrode at the bottom center is used.

Here, the molecules convect and diffuse near the surface, the rate of immobilization is determined by several factors including the inlet concentration (c0), diffusion coefficient (D), maximum surface molecule density (Γs). The convection-diffusion equation and transport-adsorption equation are used along with time-dependent study:

$${{\partial c} \over {\partial t}} + \nabla \cdot (-D\nabla c) + \textbf{u} \cdot \nabla c = R$$

$${{\partial c_s} \over {\partial t}} + \nabla \cdot (-D\nabla c_s) = k_{ads} c(\Gamma_s – c_s)-k_{des} c_s$$

For the transport-adsorption equation, it is assumed that the change of surface concentration plus the rate of surface diffusion equals the rate of Langmuir adsorption isotherm.

Fig. 6 shows the time-dependent surface concentration (cs) change when c0 = 1μM, and Fig. 7 shows the binding curve (probe density vs time) for different values of c0.

Figure 6. Time-dependent surface concentration (cs) change. t = 0~18hr.

Figure 7. Probe surface density (molecules/cm2) vs time (hr).

At concentrations above 0.1μM, the probe density almost saturates to a value of 9.6×1012 molecules/cm2 before immobilizing for 10 hours. The result highly resembles a typical binding curve, suggesting the possibility for computer simulating assisted optimization of in vitro parameters, which is really helpful for understanding underlying mechanisms for the system.

Electrochemical Impedance Spectroscopy (EIS) Simulation.

EIS is a rapid and label-free method for detection of bio-molecules, and is widely implemented on a variety of biosensors. In this task, I simulated EIS diagrams by changing the values of the heterogeneous rate constant (k0) and the double layer capacitance (Cdl). Both Cdl and k0 are affected by the immobilized surface molecule density on an electrode surface, and are important physical properties when analyzing EIS data. A slice on the xy plane is used as the modeling geometry of this task (Fig. 8).

Figure 8. Geometry being simulated for task 3. A 2D plane is sliced in the xy direction.

Here, a sinusoidal voltage wave is applied between the two electrodes (amplitude ≅ 5mV), and Bode plots and Nyquist plots can be drawn according to the measured impedance. According to the surface redox reaction, an equivalent circuit can be constructed. The equivalent circuit for this task is shown in Fig. 9.

Figure 9. Equivalent circuit used in my research and task 3.

Fick’s 2nd law and Butler-Volmer equation are used for simulation:

$${{\partial c}\over{\partial t}} = \nabla \cdot (D\nabla c)$$

$$j = nFk_0 (c_{Red} e^{ {(n-\alpha_c)F\eta} \over {RT} } – c_{Ox} e^{ {-\alpha_c F\eta} \over {RT} })$$

EIS plots are simulated for different values of k0 and Cdl (Fig. 10).

Figure 10. Bode and Nyquist plots for the simulated EIS data. k0 has a range from 0.001 to 0.1 (cm/s), and Cdl has a range from 0.01 to 100 (uF/cm2).

By undergoing this simulation project, I furthermore understood some fundamental interactions between the physical properties and outcome of my research to a new depth, and developed new concepts about how to improve it.

After completing this project, I also used COMSOL for simulating time-dependent concentration gradient variation of redox molecules in my 1st author journal paper “Diffusion impedance modeling for interdigitated array electrodes by conformal mapping and cylindrical finite length approximation”.

[Full report for this project]

Automated Microfluidic Controlling Platform

I have always wanted to automate the cumbersome experimental operation process. For biochemical experiments, often a whole day in web lab is spent in order just to acquire one single set of data. Having the hands-on experience of developing hardware and software integrated systems for several projects (e.g. Surface Plasmon Resonance Platform, Real-time Impedance Detection Systems), I initiated this project for constructing a microfluidic controlling platform that can automatically manipulate liquid-based solutions of little volume, and also assist real-time detection experiments.

Platform Structure

The platform is constructed by four sections: A Raspberry Pi, an actuator module that consists of an Arduino and two H-bridge circuits, a fluid controlling system made up of a syringe pump/syringe device, and the microfluidic platform (Fig. 1).

Figure 1. System architecture of the automated microfluidic controlling platform.

Here, the Raspberry Pi serves as the main processing unit, which an Apache web server is constructed on, and is used to communicate with a remote user by website. The front-end of the website is designed using HTML, Javascript and CSS, and the back-end is designed using PHP. Javascript and PHP communicate using jQuery, and the PHP code is written for controlling peripheral devices.

For x/y position control, the Raspberry Pi sends data through type-B USB to the Arduino, which afterwards commands the H-bridge circuits, then control the x/y stepper motors on the microfluidic platform for moving the position of the racket. The servo motor is used to control the high/low position of the tube it holds, where the high (low) position means the tube isn’t (is) inserted into the microtube (Fig. 2).

Figure 2. The microfluidic platform and surrounding modules.

For fluid control, the Raspberry Pi sends an infuse/withdraw signal to the syringe pump using another type-B USB, which subsequently pushes/pulls the syringe on it. A simple flow for moving liquid from microtube A to microtube B is:

1. Confirm that the tube (for fluid conveyance) position is high. If not, then move it up by servo motor.
2. Move to microtube A by stepper motor.
3. Change tube position to low by servo motor.
4. Withdraw liquid by syringe pump.
5. Change tube position to high by servo motor.
6. Move to microtube B by stepper motor.
7. Change tube position to low by servo motor.
8. Infuse liquid by syring pump.

Hardware Development

Fig. 3 shows a photograph of the whole platform. For the microcontrollers, Raspberry Pi 3 B+ and Arduino Uno are used. For the H-bridge, L298N dual driver module is used. Legato® 111 syringe pump (kd Scientific) and Series 700 Microliter syringe (Hamilton) are used for fluid control. HMS-25BY46L38 stepper motors and an SG90 servo motor are used for the microfluidic platform.

Figure 3. Photograph of the automated microfluidic controlling platform.

The circuitry for this platform is relatively easy compared with other projects (e.g. Aroma Alarm Clock), and is consisted simply with wires and resistors. Except the motors and the racket with microtubes, all the other components of the for microfluidic platform are fabricated using 3D printing. 3D-printed gears are fixed with the stepper motors. Precise control of racket position (~ 0.2mm) is realized by combining the gears with a 3D-printed linear gear that fits on the racket and another linear gear of a subsidiary platform which the racket sits on. Below is a clip demonstrating how the 3D-printed components, the stepper motors, the racket with microtubes, and a microfluidic electrode chip are integrated together, along with x/y position control of the racket.

Software Development

The program structure written inside Raspberry Pi is quite similar to the program structure in another project “Real-time Impedance Detection Systems”. The difference for this project is that PHP is used to directly communicate with peripheral devices (Fig. 1).

Fig. 4 displays the website-based user interface for controlling this platform. The UI is divided in four sections: procedure window, buttons field, racket window, and status window. Briefly, the user can save/load settings from the connected Arduino, add/delete a control command, start/pause/stop the current procedure, download the control procedure to a text file, and append a command after another one.

Figure 4. Website user interface of the platform.

Fig. 5 shows the settings menu. Here, the user needs to find the device url for Arduino and the syringe pump beforehand and insert them. Several controlling preferences, such as motor operation delay time, steps for the stepper motor to move per cell (microtube), precise high/low position of the servo motor, current cell of the racket … etc. can be set.

Figure 5. Settings menu of the UI.

Fig. 6 shows the add command menu. The user can either move the stepper motor to a target cell (microtube), infuse/withdraw fluid at a self-defined rate and target volume, move servo high/low position, or perform a time delay.

Figure 6. The add command menu of the UI.

Concentration Gradient Generation

An automated concentration gradient generation process is written, and is carried out by the automated platform. Here is a clip for demonstration (speed = 10x):

Real-time Impedimetric Detection

The tube doesn’t have to directly connect to a syringe. A detection chip can be inserted between for real-time detection of different solutions (similar to the method illustrated in Fig. 1 of another project Surface Plasmon Resonance Platform). A microfluidic interdigitated electrode chip using microfabrication technique is previously developed (which is a part of my research), and is used for detection of electrochemical impedimetric properties of the fluid. Here, potassium ferricyanide (K3Fe(CN)6) and potassium ferrocyanide (K4Fe(CN)6) are serial diluted using the platform, and real-time impedance detection is carried out using an electrochemical analyzer (CHI614b, CH Instruments) and the microfluidic chip.

Fig. 7 shows the detection result. It can be seen that the solution switching time is relatively fast and stable, and the detection time is consistent, which demonstrates the advantages of this automated platform.

Figure 7. Real-time impedimetric detection plot for different diluted concentrations of K3Fe(CN)6/K4Fe(CN)6.


In summary, a website-controlled automatic microfluidic controlling platform is designed and fabricated for real-time microfluidic sensing and other applications. Solution manipulation using this platform is stable, repeatable, and time-saving compared with manual operation.

Miniaturized ELISA Platform

Studying in a cross-disciplinary department meant having the freedom to choose what to explore. I entered the intelligence bio-sensing lab (previously named bio-molecular device lab) hosted by professor Lin-Chi Chen when I was a junior. There I was trained how to put into practice the engineering skills that I have learned during college, and implementing them on biosensing. I became interested in manufacturing devices that can realize automation, assist research, or help reduce the cost for lab experiments (e.g. Real-time Impedance Detection Systems, Surface Plasmon Resonance Platform, Automated Microfluidic Controlling Platform). This miniaturized ELISA platform serves as the first one among those devices and systems I had created.

The enzyme-linked immunosorbent assay (ELISA) is a commonly used analytical biochemical assay that uses antibodies against the protein to be tested to detect the presence of ligands (usually proteins) in the liquid sample. However, the traditional method for performing this assay is costly and time-consuming. Therefore, I decided to construct a miniaturized ELISA platform that can help reduce sample usage, and thus make it cheaper.

Small circular holes are cut on a thin acrylic board are by laser cut, a holder for assisting supporting the microwell is fabricated using 3D printing, and PVDF films are used as the base material for protein immobilization (Fig. 1).

Figure 1. Materials used for the miniaturized ELISA platform

The acrylic board with holes and another board with no holes are used to clip the PVDF film tight, wrapped with tape, making microwells with a volume capacity of ~10μL (Fig. 2 left). The microwell is put on the 3D-printed holder, and the right picture of Fig. 2 shows the microwell platform with each well containing 10μL deionized water.

Figure 2. PVDF clipped with an acrylic board with holes and another board without holes (left), and the microwell platform with every well containing 10μL deionized water.

Streptavidin-HRP is diluted using PBS buffer, and 5μL of the solution is added in each of the microwell. Then 5μL TMB is added for validation of the colorimetric detection method. Fig. 3 shows the experiment result using the platform for qualitative analyzing different concentrations of streptavidin-HRP. It can be seen that different concentrations yield different intensities of absorbed light signals (λ = 450nm), thus this platform can be further improved for real experimental use.

Figure 3. Different concentrations of streptavidin-HRP with TMB for colorimetric detection using the miniaturized ELISA platform.

This project is the first one for me to implement simple skills that I have learned during the first three years in university on real bio-detection research issues, which motivated me to start thinking of practical methods to use engineering techniques for solving problems in an interdisciplinary way.

Structural Biology Simulation

The usage of proteins is almost inevitable in most biochemical experiments. The ironic thing is, even if several billion or trillion proteins are present right in front of us, we never really get to see their true form due to their microscopic sizes. Thus, I enrolled in a class named structural biology, which I learned four programs: PyMOL, Swiss PDB Viewer, MolMol, and Chimera, for visualizing proteins, their physical properties, and several interaction mechanisms. This helped me understand important structural properties about the protein I had been studying.

Here I demonstrate some simulation methods implemented on the protein: signal transducer and activator of transcription 3 (STAT3). Three PDB files are used in this project: 3cwg, 1bg1, and 1bf5.

Analysis 1: Protein visualization using cartoon (top-left), dots (top-right), sticks (bottom-left) and spheres (bottom-right). Secondary structures such as the alpha helix and beta sheet are colored differently (PDB ID: 3cwg) (PyMOL).

Analysis 2: The volume of the protein (PDB ID: 1bg1) is calculated as 71.613nm3, and its surface area is calculated as 263.23nm2. The structure is transformed into a spherical molecular representation prior to calculation (Chimera).

Analysis 3: Total width and height of the protein (PDB ID: 3cwg) (Swiss PDB Viewer).

Analysis 4: Morphing between two different PDB files of the same protein (PDB ID: 3cwg and 1bg1) (Chimera). The blue structure is 3cwg, and the gray-white structure is 1bg1 in the lower figure.

Analysis 5: Electric charge on alpha helix (PDB ID: 3cwg) (PyMol).

Analysis 6: Mutation of Proline to Histidine at residue 255 (PDB ID: 3cwg) (Swiss PDB Viewer).

Analysis 7: Twisting of the ϕ and ψ angle (respectively left and right figure) at residue 255 (Proline) (PDB ID: 3cwg) (Swiss PDB Viewer).

Analysis 8: Ramachandran plots of the same protein with two different PDB files (PDB ID: 3cwg (left figure) and 1bg1 (right figure)) (MolMol).

Analysis 9: Coulomb force on protein surface. The surface is colored from red (-10kcal/mol×e) to blue (10kcal/mol×e) gradient in order to indicate differences in Coulombic forces (PDB ID: 3cwg) (Chimera).

Analysis 10: The hydrogen bond between the two SH2 domains of the STAT3 dimer (PDB ID: 1bg1) (Chimera).

Analysis 11: Morphing between STAT3 (1bg1) and STAT1 (1bf5) (another similar protein of the STAT family) (Chimera). The blue structure is STAT1, and the white structure is STAT3 in the lower figure.

Fermentation Batch Reactor

For most of what we experience in everyday life, it is rare that one can directly link the obvious outcomes with their underlying theoretical grounds. Equations and plots seem such a long distance toward their practical applications. I regard this project as an important one which links observations of a simple experiment to the complex differential equations in reaction mechanics. This mini-project comes from a homework in reaction engineering, a course I had enrolled in during college. The experiment is simple that any person can carry out using easily accessible materials. The main objective is to construct a batch reactor that can exhibit fermentation with yeast, then quantify the reactions using what we have learned on class. (~age 21, 2017)

Two commercially available sugar-sweetened beverages, glucose solution, and water are used to explored how the sugar content in them affects the fermentation rate of rapid yeast. The anaerobic fermentation of yeast in anaerobic environment is:

C6H12O6 (monosaccharide) → 2C2H5OH (ethanol) + 2CO2 (carbon dioxide) + 2ATP

In this experiment, glass containers are filled with the solutions, then instant yeast is added the each container for production of carbon dioxide. A balloon is used for trapping the gases and is used as a volume sensor, where its dimensions are measured for calculating the volume of generated CO2. The molar concentration of CO2 is calculated using the ideal gas equation PV = nRT, and the ethanol production rate is calculated by relating with the proposed reaction and using finite different method.

Figure 1. Snapshots of balloon-sealed containers with added yeast at different times.

Figure 2. Volume of balloon (Vballoon) vs time (min).

Assume an inner air pressure of P = 1atm, a temperature of T = 310K (37°C). From the ideal gas equation, the relationship between the number of moles of CO2 and its volume is n = 3.931×10-5V, which according to the reaction, also equals the number of moles of ethanol. The molar concentration of ethanol is calculated by dividing the number of moles by its volume. And by using finite difference method of the first derivative, the rate of increase for molar concentration of ethanol (rC2H5OH) is calculated (Fig. 3).

Figure 3. Increase rate of molar concentration of ethanol (rC2H5OH) vs time (min).

It can be seen that in addition to pure water (Negative), the other three sugar-containing solutions have a maximum formation rate at the beginning (marked by the blue arrow). Wherein the ethanol production rate of glucose solution is eventually lower than 0(mM/min), it is presumed either this is caused by measurement errors or that carbon dioxide is dissolved back into the liquid, causing a decrease in volume, not a decrease in the amount of ethanol.

Here the production rate of ethanol in glucose solution started at a very high value (8.29mM/min), followed by fruit tea (4.17mM/min), and then raspberry juice (3.36mM/min). However, the sugar concentration of raspberry juice is higher than that of fruit tea. There are two factors that may be affected: the type of sugar and the pH value. Among them, the pH of fruit tea is between 5.0 and 6.0 and the pH of raspberry juice is between 2.3 and 2.52. However, the optimal living environment pH of yeast is 4.5 to 5.0, so it is speculated that the acidic environment of raspberry juice inhibits the activity of yeast and reduces rC2H5OH. In addition, only glucose exists in the glucose solution, but there is sucrose in both raspberry juice and fruit tea. Sucrose can be broken down by the yeast and producing ethanol twice as much as the same concentration of glucose. This explains why the final balloon volume (408.69cm3) of fruit tea is greater than the final balloon volume of the glucose solution (361.03cm3).

As being a simple hands-on experiment, this project successfully delivered the knowledge and allowed me to learn the fundamentals through practice, by which creating a connection between reality and theory.