A python code to calculate the Brownian motion of a colloidal particle in a time varying force field. The code was developed for the research work of the Hesselink research group at Stanford University. The particle trajectory computed by the code can be useful for designing lab-on-a-chip devices. Please cite this repository and the two papers listed in the reference section when using this code.
The current example code shows the transport of a micro-particle along a moving micro-electrode array producing moving dielectrophoretic force. As the micro-electrodes are excited in sequence, the micro-particle follows the position of the active electrode.
The Brownian motion code solves the Langevin equation in discrete time. The code is general. Although the example considers a dielectrophoretic force-field, the code would work for any type of external force-field. The code would also work for nano-particles instead of micro-particles. The Brownian vibrations are more prominent for nano-particles.
The python script was tested with spyder IDE.
The Brownian motion of a colloidal particle in a low Reynolds number environment can be modeled by the Langevin equation:
$m \frac{\partial\mathbf{v}(\mathbf{r},t)}{\partial t} = \frac{kB T}{\overset{\leftrightarrow}{\mathbf{D}}(\mathbf{r})} \Bigg[\mathbf{v}_f(\mathbf{r},t)-\mathbf{v}(\mathbf{r},t) + \sqrt{2} \overset{\leftrightarrow}{\mathbf{D} {1/2}}(\mathbf{r}) \mathbf{W}(t) \Bigg] + \mathbf{F}_e(\mathbf{r},t)$
Here, $\mathbf{r} = (x_o,y_o,z_o)$ is the position of the center of the particle, $\mathbf{v}$ is the particle velocity, $\mathbf{v}_f$ is the fluid velocity, $\mathbf{F}_e$ is the external trapping/manipulation force acting on the particle, $k_B$ is the Boltzmann constant, $T$ is the temperature, $\overset{\leftrightarrow}{\mathbf{D}}$ is the diffusion tensor, and $\mathbf{W}(t)$ is a vector white noise term. The tensor $\overset{\leftrightarrow}{\mathbf{D}}_1/_2$ is defined as the element-wise square root of $\overset{\leftrightarrow}{\mathbf{D}}$. Each Cartesian component of $\mathbf{W}(t)$ is a random process unit zero mean and unit variance.
The Langevin equation can be converted into the following discrete form:
$m \frac{\mathbf{v}{i+1} - \mathbf{v}_i}{\Delta t} = \frac{k_B T}{\overset{\leftrightarrow}{\mathbf{D}}(\mathbf{r}_i)} \left[\mathbf{v} {f,i}- \mathbf{v} {i+1} + \sqrt{\frac{2}{\Delta t}} \overset{\leftrightarrow}{\mathbf{D}} {1/2} (\mathbf{r}i) \mathbf{w}_i \right] + \mathbf{F} {e, i}(\mathbf{r}_i)$
The equaiton is numerically solved using the Euler-Maruyama method.
The blue horizonal lines represent electrodes which are excited in sequence. Only one electrode is active at a time (the rest are grounded). A solid surface is assumed to be located along z = 0 plane. The red spheres are polystyrene micro-particles, each having a radius of 10 micron.
The position and velocity profiles are:
A optical trap that is active during the time interval 1s < t < 8s is simulated. The optical spot is assumed to create a Gaussian potential well. The resulting force profile is given by the gradient of the potential well. A solid surface is assumed to be located along z = 0 plane (representing the substrate). Particles hitting the solid surface experience elastic collision.
It can be noted that particles outside capturing range of the trap (i.e. beyond the range where the trapping force is significant) do not get trapped.
The position and velocity profiles are:
A large optical spot and a strong gradient force model is used here for illustration purposes. It is possible to simulate optical tweezers with arbitrary spot size and traping potential depth. A few more animated examples can be found in the Hesselink lab YouTube channel :
A microfluidic device that can sort/separate live and dead yeast cells is simulated. The device uses electrodes to apply dielectrophoretic forces. Since the material properties and hence the Clausius-Mossotti factor of live cells (read circles) differ from that of dead cells (blue circles), they experience different dielectrophoretic forces. Using this principle, the cell sorting/separation operation is accomplished.
The cell separation mechanism can be analyzed by observing the y trajectory of different cells.
This work is partially supported by the National Institute of Health (NIH) Grant R01GM138716.
- Moved the force related functions to a separate file- Interpolation function for the experimental force data has been modified accordingly- Executation/elapsed time added- Status messages added- uploaded as v1.5
- Added particle-particle elastic collision dynaimics- Works fine for two simultaneous collisions. For larger number of simultaneous collisions, some modifications are needed. -
March 2, 2021
- Decoupled the simulation time step to animation frame rate. They can be independently defined now.
March 3, 2021
- Removed excess parameters from the force functions- Made the different force functions, i.e. forceDEP.py and force_spring_trap.py more consistent with each otherOnly the import call and final time needs to be adjusted for when switching the force functions. The main outputfrom both imports are function_profile(r,t).
- Added zorder for animating the beads (beads append). This ensures that the beads are always in the foreground (compared to source geometry)- Added draw_source(t) function inside the animate() function. Also, made the draw_source() as a funciton of time for dynamicmanipulation of the source geometry (e.g. turning ON/OFF optical excitation for optical trapping demo.)
- Added x,y,z position vs time plots- Added animation for yz, and zx planes along with xy plane
- Added wall collision mechanics- Added substrate graphics- Renamed some of the modules
- Streamlined the plots- Fixed random distribution of initial position
- Steamlined the draw_geo() function within the force files- Moved the definition of ro and tfinal from the main file to the force files
- Generalized the code for working with non-homogeneous particles (particles with different radius and different mass)- Released as version 1.7.0
- Updated the Diffusion tensor + hydrodynamic interactions- Redefined ro as a (Np,) vector instead of (Np,1) vector- Updated plot font size- Released as version 1.7.2
- Bug fix: matplotlib version issues: Replaced py.sca() to py.gcf().sca()
- Added feature to plot specific frames of the animation separately (useful for saving images to file)
- Added parameters defining initial particle position range (xi_lim, yi_lim and zi_lim) inside the force functions- Extra text feature on the animation (text_string1 and text_string2)
- Improved axes limit parameter definitions
- Added fluid-velocity term
- Added particle_color() function
- Generalized wall collision dynamics
- Generalized node based geometry definition- Integration of geometry drawing with geometry definition- Integration of collision dynamics with geometry definition- Reducing the number of particle data plots to reduce clutter
- Functionalized fluid velocity
- Nt definition bug. Nt = frame_rate*40 was incorrect. It should also be multiplied by tfinal- sca() function may cause issues. It might be good be replace it with old versionexplicitly defined within geometry_draw module
Zaman, M. A., et al. “Modeling Brownian Microparticle Trajectories in Lab-on-a-Chip Devices with Time Varying Dielectrophoretic or Optical Forces.” Micromachines 12.10 (2021): 1265.
https://doi.org/10.3390/mi12101265
Zaman, M. A., et al. “Microparticle transport along a planar electrode array using moving dielectrophoresis.” Journal of Applied Physics 130.3 (2021): 034902.
https://doi.org/10.1063/5.0049126
Zaman, M. A., et al. “Controlled transport of individual microparticles using dielectrophoresis.” Langmuir 39.1 (2022): 101-110.
https://doi.org/10.1021/acs.langmuir.2c02235
This work was partially supported by the National Institute of Health (NIH) Grant R01GM138716 and 5R21HG009758.
