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]

Survival Rate Prediction Model for Startup Companies

This is an end semester project that I have done with my groupmates in college. In this study, we proposed a model to predict a startup company’s future condition using a deep multilayer perceptron (MLP) and decision tree. I mainly played the role for data quantification and literature review of this study. In our study, we built a prediction analysis model for startup companies. Furthermore, we identified what the key factors are and what influences will it have on the results. Our main hypothesis is that money, people and active days are the key factors. We built the prediction model from CrunchBase, which is the largest public database with relevant profiles about companies.

Data Quantification

For the data obtained, there are 4 groups of data that are quantified: country state, employee range, roles and countries. Regarding the states group (which only exists in USA and Canada), we think that different states have different impacts on the survival rate of a startup. A lookup table for scoring of states is defined [1][2][3]. The original data for employee quantities are a range between two numbers (e.g. 101-250). These are transformed into the average of the upper and lower bound. Employees of 10000+ are defined to be 15000. For the roles group, the data “company” is arbitrarily defined as 0.1 and “company, investor” is defined as 0.9. This is because we regard a company as wealthier and influential when it also plays a role of an investor, compared with only being a company. The other roles are transformed to 0.5. The impact of country on startup environment is also studied [4] and scores from 0.329 to 0.947 are given to countries that have above 300 startup companies in record.

Table 1. Company status and scores of 23 countries.

Country Closed Operating Acquired IPO  Score 

Neural Network Implementation

The ANN runs on a Windows 10 OS (i7-8700k CPU) with DDR4 2666MHz 16G RAM and Nvidia GTX 1070 Ti GPU. Keras is used to construct the network, and three networks of multilayer perceptron (MLP) with different number of layers are designed for comparison. There are 4 outputs of the network, indicating the probabilities of the final status which are: (1) Closed, (2) Operating, (3) Acquired and (4) IPO.

Figure 1. MLP Structure of the 3 neural networks.

The three networks have respectively 2, 4 and 6 hidden layers, with each network all starting with a 1024-neuron hidden layer and ending with a 32-neuron hidden layer (Fig. 1). Table 2 shows the training results.

Table 2. Mean square error (MSE) and categorical cross-entropy loss for the 3 networks.

LossActivation Network 1 Network 2 Network 3 
Cross-entropy loss ReLU0.6950.6910.690
Cross-entropy losssigmoid0.6920.6910.691
Cross-entropy losslinear0.6950.6920.690
Cross-entropy losstanh0.6910.6920.691

The results show that almost all results are close to 70% with few difference. In general, ReLU activation has a better result than sigmoid activation, while linear activation may outperform ReLU when the network is deep enough.

The neural network behaves like a black box. It is quite difficult to conclude significant insights just by looking at the trained parameters. However, regarding the importance for business, we used a decision tree to help us understand which factor is the most important and how important they are respectively.

Decision Tree Training

We tried three different depths of decision tree: 4 (Fig. 2), 6 (Fig. 3), and 10. We set the gain ratio to be the criterion, and the confidence to 0.1.

Figure 2. Decision tree of depth = 4.

Figure 3. Decision tree of depth = 6.

The results show that if a company is large enough to exceed 500 people, the close rates are low. In most cases, large companies are acquired by mergers and acquisitions. In addition to the number of employees, funding total amount also has a significant impact on whether a company eventually survives or closes.

Key Findings

Among the factors regarding the ability to survive of starting up companies, the factors such as the number of employees, funding total amount, and the active days have significant influences on the company’s survivability. On the contrary, the factors such as country, region or number of funding rounds do not have significant influences. Whether a company acts as an investor simultaneously will also have influences on whether the company will become an IPO or will be acquired.

Future Work

Our future work is expected to integrate the inspirative insight gained from our case study with the methodology of our own. Three items are listed as in the following:

  1. Construct a heterogeneous relationship network for survival rate prediction [5].
  2. Define a data path score according to HeteSim algorithm [6][7].
  3. Predict company survival rate using MLP, decision tree and other neural networks.
  4. Predict how much money a company will raise.


  1. Bill Murphy , The Start-up Hall of Shame (America’s 10 Worst States for Entrepreneurs), © 2018 Manuseto Ventures, inc.com/bill-murphy-jr/the-startup-hall-of-shame-americas-10-worst-states-for-entrepreneurs.html
  2. Bill Murphy , 10 Top States for Entrepreneurship and Innovation, © 2018 Manuseto Ventures, inc.com/bill-murphy-jr/ranking-the-10-top-states-for-entrepreneurship-and-innovation.html
  3. Enterprising States: States Innovate, © 2015 The U.S. Chamber of Commerce Foundation, www.uschamberfoundation.org/enterprisingstates/
  4. Zameena Mejia, The top 10 best countries for entrepreneurs in 2018, © 2019 CNBC LLC,
  5. Xiangxiang Zeng, You Li, Stephen C.H. Leung, Ziyu Lin, Xiangrong Liu, Investment behavior prediction in heterogeneous information network, Neurocomputing, Volume 217, 2016, Pages 125-132
  6. Sun, Y., & Han, J. (2012). Mining Heterogeneous Information Networks: Principles and Methodologies. Synthesis Lectures on Data Mining and Knowledge Discovery, 3(2), 1-159
  7. Shi, C., Kong, X., Huang, Y., Yu, P. S., & Wu, B. (2014). HeteSim: A General Framework for Relevance Measure in Heterogeneous Networks. IEEE Transactions on Knowledge and Data Engineering, 26(10), 2479-2492. [6702458].

Coffee Maker Alarm

This is a project which I and my classmates had finished in a course focusing on practical implementation for mechatronics and system design, and is a sequel of the previous Aroma Alarm Clock project. We designed an Android app-controllable alarm clock which can heat up objects such as a water-filled container, and intend to further improve it for making coffee, so that one can enjoy a fresh cup of coffee in the morning after awakening. Similar to the previous one, I worked as the engineer in our group, and designed and constructed the hardware system.

System Framework

Fig. 1 shows the overall system architecture, which several electronic components (e.g. clock module, LCD display…) are controlled by a microcontroller that can communicate with a mobile application using a Bluetooth module. The alarm clock can be controlled either using its own hardware input (buttons) or using the mobile app. Simply put, the system has same the functions as a typical alarm clock, but has two additional features:

1. Temperature controllable heater.
2. Bluetooth communication.

Figure 1. System architecture for the coffee maker alarm.


Arduino Uno, DS1302 chip, HC-06 and LCD1602 are used respectively for the microcontroller, real-time clock module, Bluetooth module and LCD display module. Being an end-semester project in class, the components are simply connected using a breadboard. All the hardware components except the heater are packed using patterned cardboard in a simply looking fashion (Fig. 2).

Figure 2. Hardware appearance of the coffee maker alarm.

A ceramic heater plate and temperature sensor are integrated for realizing temperature controlling (Fig. 3), and are connected out from the main hardware. The total cost for all the components is 1,741 NTD (≅ 60 USD).

Figure 3. Heater and temperature sensor.


App Inventor 2 is used as the platform for creating the Android app for controlling the hardware. Fig. 4 displays the designing interface for the mobile app using the platform.

Figure 4. Design interface using App Inventor 2.

The coding section for this platform is very unique and easy to get started. It implements a block-based programming method which developers can create procedure by dragging pre-defined blocks together to for performing a certain function. Developers with no programming background can use this kind of environment for creating mobile apps. Fig. 5 shows a gallery containing the complete code for the alarm clock controlling app.

Figure 5. Gallery of block codes for the alarm clock app.

Here’s a video demonstration for using this system:

After finishing this project, I acquired important skills for system design and mechatronics integration. This enabled me to create more interesting and sophisticated projects such as the Automated Microfluidic Controlling Platform, Remote Commandable Self-Driving Toy Car and Real-time Impedance Detection Systems.

Minesweeper AI

I got addicted to Minesweeper the first few days when I started playing this awesome puzzle game. Solving the game faster every few times gives me great satisfaction of self-fulfillment. However, for the “expert” board setting, where 99 mines are hidden within a 16×30 square grid, I had never obtained a score lower than 100 seconds. Curious of what the fastest solving rate is the one can achieve, I designed a Minesweeper AI program that can automatically play the Windows Minesweeper game.

C language is used for programming, and the algorithm flowchart is displayed in Fig. 1.

Figure 1. Algorithm flowchart of the Minesweeper AI program.

Here breadth-first-search (BFS) is used to search for uncovered squares on the grid, and a queue is used for saving uncovered squares to be analyzed. A game-playing optimized algorithm is written inside the program, and OpenCV is used to take a snapshot of the Minesweeper window region and process the image for subsequent calculations. Mouse movement and click actions are realized by including the windows.h header.

The usage of the program is relatively simple. The user should only execute the Minesweeper game program, and modify the game setting beforehand. After opening the AI .exe file, the program will automatically locate the Minesweeper game, calculate the dimensions, then start playing. If a mine is accidentally clicked, the program will continue to play the next game until a fully uncovered board is achieved.

Here is a clip demonstrating the Minesweeper AI playing an expert level game and winning in 7 seconds. It failed on the first try almost at the end, but succeeded on the second try.

[Download program for the Minesweeper AI (Can only run on a windows 64-bit OS)]

Bladeless Fan

The Dyson cooling fan is an eye-catching product. At first sight, some people may wonder how the seemingly bladeless fan really works because it simply looks like a structure with no air outlet. The fact is it can does have a small outlet at the inner part of its “ring”, and has the ability to take fluid dynamics into practice and enhance the air flow, making it also an air multiplier.

Figure 1. Dyson bladeless fan and simulated air flow.

In fluid dynamics class, I and a classmate of mine decided to construct a bladeless fan by our own and study the air multiplying phenomenon. We used 3D printing to fabricate the “ring” part of the fan, and a small centrifugal fan for connecting with the inlet of the ring part and inject a strong current inside.

Figure 2. Illustration and photograph of the bladeless fan.

Now we wanted to simply test whether this structure really leads to an air multiplying effect. We divided the outlet are into 9 sections, which can be represented using a 3×3 rectangular grid, and calculated the wind speed of every section at different input voltages for the centrifugal fan.

Fig. 3 shows the air velocity profile of the output wind, and Fig. 4 shows the magnification of air flow.

Figure 3. Air velocity (t_m [=] m/s) of 3×3 section grid at different input voltages.

Figure 4. Magnification of air flow at different input voltages of the centrifugal fan.

During this project, I learned more about the fundamentals of fluid mechanics, and memorized the relevant rules more deeply, which made me have a better understanding of this subject.