Passive and active colloidal chemotaxis in a microfluidic channel: mesoscopic and stochastic models

TL;DR

This paper uses mesoscopic and stochastic models to demonstrate and analyze chemotaxis of colloidal particles in microfluidic channels, providing insights into active and passive chemotactic behaviors.

Contribution

It introduces a particle-based simulation framework and a stochastic model to explain colloidal chemotaxis in microfluidic environments, extending previous simpler models.

Findings

Chemotactic forces depend on chemical concentration gradients.

Passive and active colloids exhibit distinct chemotactic responses.

The stochastic model successfully explains experimental and simulation results.

Abstract

Chemotaxis is the response of a particle to a gradient in the chemical composition of the environment. While it was originally observed for biological organisms, it is of great interest in the context of synthetic active particles such as nanomotors. Experimental demonstration of chemotaxis for chemically-powered colloidal nanomotor was reported in the literature in the context of chemo-attraction in a still fluid or in a microfluidic channel where the gradient is sustained by a specific inlet geometry. In this work, we use mesoscopic particle-based simulations of the colloid and solvent to demonstrate chemotaxis in a microfluidic channel. On the basis of this particle-based model, we evaluate the chemical concentration profiles in the presence of passive or chemically active colloids, compute the chemotactic force acting upon them and propose a stochastic model that rationalises our findings on colloidal chemotaxis. Our model is also able to explain the results of an earlier simulation work that uses a simpler geometry and to extend its interpretation.

Figures

Tables

Table 1. Table 1: Simulation parameters for the chemotactic cell. The parameter files for reproducing the simulations are available publicly, see appendix  A .

Parameter	symbol	value
Number density	$n$	10
Solvent particle mass	$m$	1
MPCD fluid density	$\rho$	10
Cell size	$a$	1
MPCD time step	$\tau$	0.5
MD steps per MPCD step	$N_{MD}$	50
Cell dimensions	$L_x,L_y,L_z$	90, 60, 15
Buffer length	$L_{\mathrm{buffer}}$	20
Temperature	$k_{B}T$	$1/3$
Acceleration	$g$	0.001
MPCD collision angle	$\Theta$	2.6
Reaction probability	$p$	1
Colloid radius	$\sigma$	3
Unit interaction energy	${ϵ}_A$	1

Table 2. Table 2: Characteristics of the microfluidic channel and of the fluid flow, both in the experiment of Ref.  [ 9 ] and in our simulations.

Number	Simulation	Experiment
Width ($L_y$)	60	600$\mu$m
Height ($L_z$)	15	30$\mu$m
Average flow velocity $v_{av}$	0.063	2.16mm/s
Maximum flow velocity $v_{max}$	0.095	3.24 mm/s
$\mathrm{Pe}$	14	32
Re	0.48	0.076

Equations

$$r_i(t+dt)=r_i(t)+v_{i}dt+\left(g+f_i(t)\frac{d{t}^2}{2m_i}\right),$$

$$v_i(t+dt)=v_i(t)+\frac{dt}{2m_i}\left(f_i(t)+f_i(t+dt)+2m_{i}g\right),$$

$$dt=\frac{\tau }{N_{\text{MD}}},$$

$$v_i^{\prime }=v_{\xi }+{\Omega }_{\xi }\left(v_i-v_{\xi }\right)$$

$$V_{ij}=\left\{\begin{array}[]{l l}4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}+\frac{1}{4}\right)&\textrm{for }r_{ij}\leq\sigma_{ij}\times 2^{1/6}\\ 0&\textrm{else}\end{array}\right.$$

$$V_i={ϵ}_{\mathrm{wall},i}\left\{\frac{3\sqrt{3}}{2}\left({\left(\frac{{\sigma }_{\mathrm{wall},i}}{z_i}\right)}^9-{\left(\frac{{\sigma }_{\mathrm{wall},i}}{z_i}\right)}^3\right)+1\right\},$$

$$A+C\to B+C$$

$$\eta =$$

$$+\frac{m}{18a\tau }(n-1+e^{-n})(1-\cos \Theta ),$$

$$D=\frac{k_{B}T\tau }{2m}\left(\frac{3n}{(n-1+e^{-n})(1-\cos \Theta )}-1\right).$$

$$v_x(z)=v_{max}\frac{z(L_z-z)}{(L_z/2{)}^2}=\frac{\rho g{L}_z^2}{8\eta }\frac{z(L_z-z)}{(L_z/2{)}^2}.$$

$$\mathrm{Pe}=L_z{v}_{av}/D,$$

$$\mathrm{Re}=v_{\text{flow}}{L}_z/\eta ,$$

$$\mathrm{Ma}=\frac{v}{v_{\mathrm{sound}}},$$

$${\partial }_x{c}_{\alpha }(x,y)=\frac{D}{v_{\text{flow}}}{\partial }_y^2{c}_{\alpha }(x,y)$$

$$c_A(x,y)$$

$$c_F(x,y)$$

$$c_B(x,y)$$

$$\sum _{\alpha }{c}_{\alpha }(x,y)=\rho ,$$

$$\left\{\begin{array}[]{l l}x&=r\cos\varphi\sin\theta\cr y&=r\cos\theta\cr z&=r\sin\varphi\sin\theta\end{array}\right.$$

$$R{k}_D{\partial }_r{c}_A=k_0{c}_A,$$

$$\vec{F}=\sum _{\alpha }\int d\mathbf{r}{c}_{\alpha }(\mathbf{r})\vec{\nabla }{V}_{\kappa ,\alpha }(\vec{r}).$$

$$\vec{F}$$

$$=\sum _{\alpha }{\int }_{r=R}d\mathbf{r}{c}_{\alpha }(R\hat{r}){\vec{1}}_r\int dr{e}^{-\beta V_{\kappa ,\alpha }(r)}{\partial }_r{V}_{\kappa ,\alpha }(r),$$

$$\vec{F}=\frac{2}{\beta }\sum _{\alpha }{\Lambda }_{\kappa ,\alpha }{\int }_{r=R}d\mathbf{r}{c}_{\alpha }(R\hat{r}){\vec{1}}_r,$$

$${\Lambda }_{\kappa ,\alpha }={\int }_0^{R}drr\left(e^{-\beta V_{\kappa ,\alpha }(r)}-1\right).$$

$${\vec{F}}_N=\frac{2}{\beta }\left({\Lambda }_{N,A}-{\Lambda }_{N,F}\right){\int }_{r=R}d\mathbf{r}{c}_A(R\hat{r}){\vec{1}}_r,$$

$$F_{N,y}$$

$$=\frac{8\pi }{3\beta }\left({\Lambda }_{N,A}-{\Lambda }_{N,F}\right)\lambda R,$$

$$R{k}_D{\partial }_r{c}_A=k_0{c}_A+R{k}_D\lambda \cos \theta$$

$$c_A=c_0+c_1\frac{R}{r}+c_2{\left(\frac{R}{r}\right)}^2\cos \theta +\lambda r\cos \theta .$$

$$\left\{\begin{array}[]{l l l}c_{0}&=c_{A}(x,y)\\ c_{1}&=-\frac{k_{0}}{k_{0}+k_{D}}c_{0}\\ c_{2}&=-\frac{k_{0}}{k_{0}+2k_{D}}\lambda R\end{array}\right.~{},$$

$$c_B=-c_1\frac{R}{r}-c_2{\left(\frac{R}{r}\right)}^2\cos \theta .$$

Full text

1 Introduction

Nanomotors are nano- to micro-meter sized machines that use a local (internal or found in the close environment) source of energy to move or perform work. The catalytic conversion of a chemical fuel on catalytically coated colloids was used with success for different types of motors (metallic rods [1], polystyrene [2] and silica [3] Janus particles, dimers [4]). Chemically powered nanomotors represent very promising devices for the execution of tasks such as sensing or cargo delivery in nano- to micro-meter scaled environments [5, 6, 7]. Understanding the response of nanomotors to chemical concentration gradients, chemotaxis, is critical to engineer these possible applications.

The experimental characterisation of chemotactic motion can proceed either via chemical sources that supply the environment with fuel, at specific “target” locations, or via a flow that allows the sustainment of the gradient. The first idea was used by Hong et al to demonstrate the occurrence of chemotaxis for rod nanomotors [8]. The second strategy was used by Baraban et al in Ref. [9] where the authors studied experimentally the chemotactic motion of chemically powered nanomotors in a microfluidic channel and observed a lateral deviation of the nanomotors when the fuel is input asymmetrically at the inlet of the cell. The successful modeling and understanding of colloidal chemotaxis is also relevant for biological systems, with the traditional example being the chemotactic behaviour of bacteria [10]. The focus on biological chemotaxis has been put forward by Sengupta et al who also used a microfluidic channel to observe the enhanced migration of enzymes across a microfluidic channel in the presence of a substrate gradient [11]. More recently, the sorting of enzymes on the basis of their chemotactic response has been demonstrated in Ref. [12].

Numerical simulations of chemotaxis can either assume an expression for a “chemotactic force” or reproduce the chemotactic process itself, using a direct simulation of the solvent. This article starts with the latter approach, using Molecular Dynamics simulation with an explicit solvent in which several chemical species are represented. We then formulate a continuous picture for the solvent concentration field and use it to build a stochastic model that captures the chemotactic behaviour and its connection to our simulation parameters.

The literature on particle-based modeling of chemotaxis for nanomotors is rather scarce. Chen et al [13] considered a simplified system in which a constant chemical gradient is imposed by the boundaries of the simulation cell. In this work, we seek to imitate the experimental setup of Baraban et al [9] and to improve the theoretical understanding of this experimentally relevant configuration for chemically powered nanomotors, albeit using a simpler dimer-type nanomotor instead of the Janus and tubular microjet motors in the experiment. Our computational experiments builds gradually by starting with a non-catalytic spherical colloid, then adding a catalytic property to the spherical colloid and finally using the dimer nanomotor. We come back to the results of Chen et al [13] and extend their conclusions.

In section 2, we review the mesoscopic simulation model. The simulational implementation of the experiment of Ref. [9] is laid out in section 3. The continuous representation of the fluid’s diffusion profile, including in the situation where chemical reactions occur on the surface of the colloid, is given in section 4. There, we also compute the chemical gradient induced forces on the colloids. The stochastic model for the colloids is presented in section 5. We give the results for both modeling strategies in section 6 and conclude in section 7. In appendix A, we provide full information on how to access the numerical code to reproduce our findings.

2 Simulation model

The solvent consists of point particles with a mass $m_{i}$, a position $r_{i}$ and velocity $v_{i}$. There is no force acting between solvent particles, their interaction is instead modeled via cell-wise collisions at fixed time intervals $\tau$ using the Multiparticle Collision Dynamics (MPCD) collision rule introduced by Malevanets and Kapral [14, 15]. MPCD has been used successfully to investigate the dynamics of colloids in microfluidic channels by Prohm et al, in the context of inertial focusing [16], or by Nikoubashman et al to study the flow of colloids in the presence of obstacles [17], for instance.

The evolution of fluid particles, in the presence of forces due to colloids or to the flow-inducing field, is resolved numerically using the velocity Verlet algorithm [18, 19]

[TABLE]

[TABLE]

where $f_{i}$ is the total of pair forces on particle $i$ and $g$ is the external acceleration. The timestep $dt$ for MD is a fraction of the one for MPCD

[TABLE]

where $N_{\textrm{MD}}$ is the number of MD steps between successive MPCD collisions.

At fixed time intervals $\tau$, the fluid particles are sorted in a lattice of cubic cells of side $a$, whose origin is shifted randomly to ensure Galilean invariance [20], and their velocities are collided cell-wise according to

[TABLE]

where the prime denotes the post-collision value, $\xi$ is the cell containing the particle $i$, $\Omega_{\xi}$ a rotation operator of angle $\Theta$ in $\mathbb{R}^{3}$ around a randomly chosen axis and $v_{\xi}$ is the centre-of-mass velocity of the cell. The transport properties of a MPCD fluid can be computed analytically [21, 22, and references therein]. The simulation parameters are given in table 1.

One technique to obtain a Poiseuille flow is to impose a constant acceleration field in the simulation cell [23], in combination with periodic boundary conditions in the direction of the flow ($x$) and stick boundary conditions in the direction transverse to the flow ($z$). In Refs. [24, 19], this technique was applied to MPCD fluids. In the $y$ direction we use specular boundary conditions for the fluid, so that the flow velocity profile $v_{x}(z)$ bears no dependence on $y$. We implement stick boundary conditions for the flow in the $z$ direction with a combination of bounce-back collisions and ghost particles at the boundary MPCD cells [25]. The ghost particles also set the temperature at the wall. In this work the walls also provide a sufficient thermostatting action to compensate for the flow-induced heating and no bulk thermostatting is used. We have verified that the temperature of the fluid remains stable after a transient period at a value that is about 3% in excess of the temperature set at the walls.

The colloids evolve according to the velocity Verlet algorithm, similarly to the solvent particles, but do not participate in the collisions and are not subject to the acceleration field $g$. In the case of dimer nanomotors, the distance between the two spheres is held constant using the RATTLE algorithm [26].

The coupling with solvent particles is done via the shifted and truncated Lennard-Jones 12-6 potential, of the form

[TABLE]

where $\epsilon_{ij}$ denotes the strength of the potential and $\sigma_{ij}$ the radius of the colloid. $i$ and $j$ represent the species of the solvent and colloidal particles and $r_{ij}$ is the distance between them. The variation of $\epsilon_{ij}$, depending on what type of solvent and colloid interact, leads to a net force in the presence of chemical concentration gradients. This will be made explicit in section 4. $\epsilon_{ij}$ for interactions with the solvent of type $A$ is set to $1$ and defines the energy scale. In the following, all quantities are expressed in simulation units with length $a$ of the unit cell, mass $m$ of the solvent particles, energy $\epsilon_{A}$ and time $\sqrt{a^{2}m/\epsilon_{A}}$. As the interaction between all colloids and solvent particles of species $A$ is equal, the first subscript to $\epsilon_{\kappa,A}$ is dropped in the following.

Confinement of the colloids in the $z$ direction is obtained by a purely repulsive Lennard-Jones 9-3 potential, of the form

[TABLE]

where $z$ is here the distance to the closest horizontal wall. The same potential is applied in the $y$ direction to avoid hypothetical crossings of the lateral wall during a simulation.

MPCD fluids are well suited to describe chemical reactions influenced by catalytic surfaces [27] or in the bulk [28]. At the surface of a catalytic colloid $C$, the reaction

[TABLE]

converts fluid particles of type $A$ (the fuel) to fluid particles of type $B$ (the product). The reaction occurs when the fluid particles crosses the interaction region of the colloid and is executed with a probability $p\in[0,1]$, when the fluid particles exits the interaction region to avoid any discontinuity in the energy [27].

All simulations were performed with the open-source RMPCDMD software [29, 30].

3 Cell design

The inspiration for the design of the microfluidic cell comes from the experimental work on the chemotaxis of nanomotors by Baraban et al [9]. There, a channel is fed with fluid at a fixed flow rate via three inlets. As only an inlet contains a chemical species of interest, i.e. the fuel for the nanomotors, the lateral distribution of species is inhomogeneous. The effect of the resulting gradient is a chemotactic behaviour that is observed by monitoring the deviation angle of the nanomotors with respect to a straight line motion.

To reproduce the experimental features in a simulation, we have setup a thin channel with a forced flow in the program chemotactic_cell of the RMPCDMD software that is illustrated in Fig. 2. A simulation snapshot in Fig. 2 shows the Poiseuille flow, the concentration field for the fluid species $A$ and the trajectory of a dimer nanomotor. To avoid reproducing the inlet channels at a large computational cost, we assign the solvent species to be either $A$ (the fuel) or $F$ (a neutral fluid species) arbitrarily, depending on the lateral position of the particles. The total number of particles in the simulation is constant and particles are not created nor destroyed but “recycled” when they re-enter the simulation cell due to the periodic nature of the $x$ boundary. The region of the cell where the forced attribution of species is applied is called the buffer and is located in the region $0<x<L_{\mathrm{buffer}}$. In the following of the paper, the trajectories from the mesoscopic simulations are shifted in $x$ by $-L_{\mathrm{buffer}}$ to match the coordinates system of the stochastic model.

Although the simulations are three dimensional, the motion of the colloids is limited around the centre of the cell in the z-direction by the confining walls.

The flow between two plates is of the Poiseuille type with maximal velocity $v_{\textrm{flow}}$ and average velocity $v_{av}=2/3~{}v_{\textrm{flow}}$. We report the characteristic numbers of the flow in Table 2. For the work of Baraban et al [9], we use the values found in the article for the flow rate (140$\mu$L per hour), the width (200$\mu$m per inlet, there are three inlets) and temperature (300K). The height of the channel was confirmed by email to be 30$\mu$m. For the properties of water, we use reference data from the NIST [31, retrieved January 11, 2017] (density $\rho_{\mathrm{water}}=996.56$ kg/m3, and viscosity $\eta_{\mathrm{water}}=8.54~{}10^{-4}$Pa$\cdot$s) and take the diffusion coefficient $D_{exp}=2~{}10^{-9}m^{2}/s$ that is a reasonable approximation for both water [32] and hydrogen peroxide [33]. For the MPCD fluid properties, we refer to the review by Kapral [22]. The viscosity is

[TABLE]

and the self-diffusion coefficient is

[TABLE]

We have verified that the velocity profile obeys the theoretical prediction [24]

[TABLE]

We have aimed for a similar fluid regime with respect to the experiment, that is a high Péclet number ($\mathrm{Pe}$) flow, as our estimate for the concentration profile relies on $\mathrm{Pe}\gg 1$ (see section 4 and Ref. [34]). A value as high as for the experiment could not be obtained, but the important feature is that the transport by the flow dominates the one by diffusion in the $x$ direction. The Reynolds number $\mathrm{Re}$ should remain moderate for the laminar regime to hold.

The Péclet, Reynolds, and Mach number for the flow are computed as

[TABLE]

[TABLE]

and

[TABLE]

where $v_{\mathrm{sound}}=\sqrt{\frac{5}{3}\frac{k_{B}T}{m}}$ is the speed of sound [16]. The maximum Mach number is obtained for the maximum velocity of the flow $v_{max}=\frac{\rho~{}g~{}L_{z}^{2}}{8\eta}\approx 0.095$ and is $\mathrm{Ma}\approx 0.13$.

Besides the geometry of the cell, the simulation protocol differs slightly from its experimental counterpart. The colloid moves on a track at fixed $y=\frac{L_{y}}{2}+y_{\mathrm{shift}}$ and $z=\frac{L_{z}}{2}$. This restriction is lifted when the $x$ position of the colloid (centre of mass) has passed $L_{\mathrm{buffer}}$ plus its own radius. This allows for a systematic comparison of the chemotactic drift across repeated runs without suffering from possible disturbance from the resetting of particles in the inlets. The shift $y_{\mathrm{shift}}$ in the $y$ direction ensures that solvent particles in the interaction range of the colloid are only of species $F$ and not influenced to a change of species in the buffer region.

4 Density profiles and surface interaction

In this section, we compute the stationary concentration field $c_{\alpha}(\mathbf{r})$ for the different chemical species in the cell and the resulting chemotactic force on the colloids.

We approximate the evolution of $c_{\alpha}(\mathbf{r})$ by a 1D diffusion equation where the spatial direction of the flow is proportional to time, i.e. $x=v_{\textrm{flow}}t$, and the diffusion process acts on the transverse direction $y$. This approximation has been tested experimentally in Ref. [34] and is only valid close to $z=L_{z}/2$. Close to the boundaries $z=0$ and $z=L_{z}$ of the cell, a different transport regime is taking place. Here, we neglect the variations in the $z$ direction for the computation of the concentration $c_{\alpha}(\mathbf{r})$ in the absence of catalytic particle.

Fixing $z=L_{z}/2$ and identifying the time and the $x$ coordinate results in

[TABLE]

The separate inlets for the $A$ and $F$ chemical species translate into the initial condition $c_{A}(0,y)=c_{0}\Theta(L_{y}/2-y)$, $c_{B}(0,y)=0$ and $c_{F}(0,y)=c_{0}\Theta(y-L_{y}/2)$ for $x=0$ (equivalently $t=0$). The solution to Eq. (12) is

[TABLE]

where $\operatorname{erfc}$ is the complementary error function. In the absence of a catalytically coated colloid, the value of $c_{B}$ remains zero at all times. The average number density is a constant,

[TABLE]

in the low Mach number conditions here.

In the remainder of this section, we use coordinates centered around the colloid. The spherical coordinates are defined by the following relations

[TABLE]

and are represented in Fig. 3.

The presence of a catalytic colloid is taken into account in Eq. (12) by a radiation boundary condition (RBC) on the surface of the colloid, at radius $R$. The boundary condition is applied at the limit of the interaction region, i.e. $R=2^{1/6}\sigma$, where the continuum diffusive picture breaks down and where the chemical reactions are triggered. The RBC identifies the flux of chemicals with the consumption of the catalytic reaction, at the surface of the colloid. It was developed by Collins and Kimball [35] and rederived later by other authors [36, 37]. The RBC at radius $R$ in the absence of external gradient, for the species $A$, is

[TABLE]

where $k_{D}=4\pi RD$ is the diffusion-limited rate constant and $k_{0}=pR^{2}\sqrt{8\pi k_{B}T/m}$ is the chemical rate [38]. $p$ is the reaction probability defined in section 2 and $m$ is the mass of a fluid particle.

4.1 Surface interaction

Following the work of Rückner and Kapral [27], we sum the interaction potential $V_{\kappa,\alpha}$ between a colloid of species $\kappa$ and the fluid particles of species $\alpha$ over the interfacial region to obtain the potential energy. Differentiation with respect to the colloid’s coordinate gives the force on the colloid:

[TABLE]

In all the subsequent derivations, we rely on the solution of the diffusion equation in the bulk of the fluid. Within the interfacial region, the interaction with the colloid modifies the concentration $c_{\alpha}$ by the Boltzmann weight [27]. Equation (19) then reduces to

[TABLE]

where $\beta=\left(k_{B}T\right)^{-1}$ and $\vec{1}_{r}$ is the unit vector pointing in the radial direction. Integrating by parts yields

[TABLE]

where we have defined

[TABLE]

Within this framework, it is sufficient to compute surface integrals over the colloids to compute the force acting on them. The use of Eq. (22) in the literature is however limited to situations where there is no external chemical gradient.

4.2 Single passive colloid

For a single passive colloid of type $N$, we use the solution (13)-(15) together with Eq. (22).

[TABLE]

Locally around the colloid, we use a first-order expansion of $c_{\alpha}$ in the lateral direction.

[TABLE]

where $\lambda=\partial_{y}c_{A}(x,y)$ The explicit dependence of $F_{N,y}$ and $\lambda$ on the coordinates is omitted for brevity.

4.3 Single active colloid

To take into account the catalytic activity of a colloid, the diffusion equation (12) is solved with the radiation boundary condition (RBC). Equation (18) expresses the flux that originates from the chemical activity of the colloid, that is the only contribution to the flux in the situation where the RBC was derived. The total radial flux is the sum of the reaction-induced flux from Eq. (18) and of the diffusive flux444The flux is multiplied by the surface of the sphere $4\pi R^{2}$ as this is how the RBC is expressed in the literature. $4\pi R^{2}D\vec{1}_{r}\cdot\vec{\nabla}c_{A}=Rk_{D}\lambda\cos\theta$.

[TABLE]

We make the following ansatz for the concentration field $c_{A}$:

[TABLE]

Eq. (28) is a solution of the diffusion equation and matches the boundary condition (27). The first three terms come from a truncation of the Legendre polynomial expansion of a diffusion profile with a spherical catalytic sink. The last term is needed to reflect the presence the external gradient. We obtain the coefficients by inserting Eq. (28) in Eq. (27):

[TABLE]

where we have used for $c_{0}$ the solution (13). The value of $\lambda$ is obtained by differentiating Eq. (13) with respect to $y$.

It is important to mention that (i) in the absence of a gradient, the coefficients in Eq. (29) lead to the standard solution of a spherical sink and (ii) in the absence of chemical activity, the solution (13) is recovered to linear order.

To obtain the solution for $c_{B}$, we observe that the reactive flux of $A$ particles absorbed at the surface of the sphere equals the one of $B$ particles, except that the currents flow in opposite directions. The solution to the diffusion equation with these opposite flows on the surface is

[TABLE]

To compute the force in Eq. (22), we set $\Lambda_{C,F}=\Lambda_{C,A}$. Using Eq. (16), we find that the chemotactic force on the active colloid $C$ is

[TABLE]

4.4 Dimer nanomotor

A dimer nanomotor is made of two spheres linked rigidly, with a distance $d$ between their centre of masses. One ($C$) is catalytic and acts as a sink for $A$-species solvent particles and a source of $B$-species solvent particles, creating locally a gradient centered around $C$.

For the dimer nanomotor, the concentration fields depend on the presence of the catalytic sphere of the motor in the same way as for the single active sphere and we reuse Eqs. (29) and (31). The force must be evaluated also on the non-catalytic sphere $N$ where the spherical symmetry does not hold. We evaluate this expression numerically as

[TABLE]

where the prime denotes the spherical coordinates around the $N$ sphere of the dimer, illustrated in Fig. 4.

The forces on $C$ and $N$ are summed to obtain the centre-of-mass ($\mathrm{com}$) force:

[TABLE]

The torque on the dimer is

[TABLE]

The motion of the dimer is studied via its centre-of-mass position in the $x$-$y$ plane and its inclination $\phi$ with respect to the $x$ axis.

5 Langevin dynamics for the colloids

For the spheres, the evolution of $x$ and $y$ (here, coordinates in the laboratory frame of reference) is given by the overdamped Langevin equations [39].

[TABLE]

where $v_{\textrm{flow}}$ is the average velocity due to the flow at $z=L_{z}/2$, $\xi_{x}$ and $\xi_{y}$ are normal distributed white noise, $D$ is the diffusion coefficient of the colloid and $\gamma$ its friction. For the single spheres, we use $\gamma=4\pi\eta\sigma$, the slip hydrodynamics friction, and $D=k_{B}T/\gamma$.

For the dimer, we use the following Langevin equations for the centre of mass and orientation

[TABLE]

where the projected forces $F_{\parallel}$ and $F_{\perp}$ are

[TABLE]

and $\mathcal{T}$ is the torque defined in Eq. (34). The rotation operations in Eq. (44) and (46) are due to the difference in friction parallel and transverse to the axis of the dimer.

The diffusion coefficients for the dimers are obtained by performing equilibrium simulations in a cell of dimension $\vec{L}=(32,32,15)$ and with no chemical reaction, all other simulation parameters being taken equal. The cartesian and angular velocity distributions for the dimer are shown in Fig. 5.

The dimer is an anisotropic colloid with axial symmetry. Accordingly, we compute separately the diffusion coefficient parallel ($\parallel$) and transverse ($\perp$) to its axis as

[TABLE]

where $\zeta$ is $\parallel$ or $\perp$, following Ref. [39]. The projected velocities are defined by

[TABLE]

with $\hat{u}$ the unit vector in the direction joining the two spheres. The angle $\phi$ measures the deviation of $\hat{u}$ in the x-y plane with respect to the unit vector $\vec{1}_{x}$. The results of these simulations is $D_{\parallel}=2.0~{}10^{-3}$ and $D_{\perp}=1.5~{}10^{-3}$. The rotational motion is characterised in the same fashion, using the velocity autocorrelation of the angle $\phi$: $\langle\dot{\phi}(\tau)\dot{\phi}(0)\rangle$ with the result $D_{r}=1.4~{}10^{-4}$. We determine the friction coefficient using the fluctuation-dissipation relation

[TABLE]

The rotational time $D_{r}^{-1}\approx 7100$ is such that in the chemotaxis simulations the dominant angular behaviour reflects the rotational drift. The contribution of rotational diffusion is nevertheless visible in the spread of angles in Fig. 10.

6 Results

6.1 Passive sphere

Here, we place a single sphere at the entry of the cell, in the upper y inlet surrounded by solvent particles of species $F$. When $\epsilon_{A}\neq\epsilon_{N,F}$, the y component of the force due to the chemical gradient is non-zero.

A typical trajectory for the passive sphere is represented in the x-y plane in Fig. 6. The main component is a displacement to the right under the influence of the flow, to which thermal fluctuations are superimposed. For this trajectory, $\epsilon_{N,F}=\epsilon_{A}$ so there is no chemotactic behaviour.

Upon changing $\epsilon_{N,F}$, a lateral force is exerted on the colloid due to the combination of the asymmetry of $c_{A}$ and $c_{F}$ with respect to the colloid and of the difference in surface interaction between the colloid and the $A$ and $F$ solvent species. This situation is one of passive diffusiophoresis.

The results for passive spheres are summarised in Fig. 8 for all chosen values of $\epsilon_{N,F}$. For the central panel $\epsilon_{N,F}=\epsilon_{A}=1$, the sphere moves in $x$ down the flow (i.e. from left to right in the figure) and undergoes diffusive motion in $y$. For $\epsilon_{N,F}<\epsilon_{A}$, we have $\Lambda_{N,A}-\Lambda_{N,F}<0$ and we expect from Eq. (26) a lateral force whose sign is opposite to the gradient $\lambda$. As $c_{A}(x,y)$ decreases for increasing $y$, $\lambda<0$ and the chemotactic force is upward, as confirmed by the panels $\epsilon_{N,F}=0.25$ and $\epsilon_{N,F}=0.50$. Conversely, for $\epsilon_{N,F}>\epsilon_{A}$, the chemotactic drift is negative. This is shown for $\epsilon_{N,F}=2$ and $\epsilon_{N,F}=4$ in Fig. 8.

Figure 8 reports the results for both the mesoscopic model and the stochastic model. Given the fluctuating nature of the colloid dynamics, we have repeated the simulations 16 times and denote by a shaded area one standard deviation below and above the average results. We observe that the average trajectories are similar and also that the spread of trajectories is similar.

This first result on passive chemotaxis provides here the simplest context for the experimental setup and allows a possible calibration of the quantity $\Lambda_{N,A}-\Lambda_{N,F}$ without any dependence on reaction rates.

6.2 Active sphere

We now turn to the chemotactic behaviour of an active sphere. The injection setup is the same as for the the passive sphere but now $\epsilon_{C,F}=\epsilon_{A}$ and the mechanism that creates a systematic lateral motion for the passive sphere, that is proportional to $\Lambda_{C,A}-\Lambda_{C,F}$, is effectively zero.

What occurs, instead, is that the local gradient in $c_{A}$ generates an asymmetric distribution also for $c_{B}$ that, together with a nonzero value for $\Lambda_{C,B}-\Lambda_{C,A}$, creates an original combination of passive and active diffusiophoresis. Using Eq. (31) and observing that $c_{2}>0$ (as $\lambda<0$ everywhere in the simulation), we expect a downward chemotaxis for $\Lambda_{C,B}-\Lambda_{C,A}>0$ (i.e. $\epsilon_{C,B}<\epsilon_{A}$) and an upward chemotaxis for $\Lambda_{C,B}-\Lambda_{C,A}<0$ (i.e. $\epsilon_{C,B}>\epsilon_{A}$). We explore the effect of changing $\epsilon_{C,B}$ in Fig. 8 that confirms this direction for the chemotaxis of the active colloid. As for the passive sphere, the comparison between the mesoscopic and the stochastic models is positive.

Even though the situation of the active sphere is simpler than the experiments on nanomotors of Ref. [9], it already contains a complex ingredient: the chemotactic force is entirely caused by the self-generated concentration field around the colloid. Thanks to the derivations in section 4, we understand how the asymmetry of the imposed concentration fields of $A$ and $F$ lead to an asymmetry in the concentration of $B$ that gives rise to chemotaxis.

6.3 Dimer nanomotor

For the dimer nanomotor, we track in the simulations the centre-of-mass position $\vec{r}_{\mathrm{com}}$ and the orientation $\phi$. The release from the injection track occurs here as soon as both spheres have exited the buffer region and the initial orientation is $\phi=0$.

We have chosen, for the simulation parameters, that $\epsilon_{C,\alpha}=\epsilon_{N,\alpha}$, for all solvent species $\alpha$. This is at variance with simulations of the prototypical dimer nanomotor [27, 40, 4, 13] where $\epsilon_{C,A}=\epsilon_{C,B}=\epsilon_{A}\neq\epsilon_{N,B}$. Preliminary tests showed that in this situation the trajectories are pathological in the sense that there was no “gentle” deviation of the orientation $\phi$ and that no systematic tendency could be found. For nanomotors in a bulk environment with no a priori gradient in the chemical concentrations, this choice would be irrelevant.

Both the lateral deviation in the x-y trajectories and the angular deviation for $\phi$ are reported in Figs. 10 and 10 and are consistent between the mesoscopic simulations and the stochastic model. This confirms the adequateness of the combined force-torque evolution given by Eqs. (44)-(45) as both the direction in $y$ and in $\phi$ are reproduced. The origin of the lateral deviation can be understood as coming not only from the redirection of the self-propelling force, that is oriented along the dimer axis, but also from the net lateral force on the centre of mass whose direction we can infer from Eqs (29), (31) and (33). To our knowledge, this is the first time that the continuum force computation of Ref. [27] is extended to compute the lateral force on a nanomotor and the resulting torque on the dimer nanomotor.

6.4 Comparison to earlier simulation work

Aside from the geometric considerations that are specific to the microchannel, the stochastic model that we have designed can accommodate other situations. We use it to test, in a qualitative sense, the observations of Ref. [13]: we re-use the diffusion coefficients that we have obtained with our cell geometry and fluid parameters but change the type of gradient and the interaction type of the motor to match those of the earlier paper, which we call the “constant gradient” model. In Ref. [13], Chen et al design a cell in which the gradient of the chemical species $A$ and $F$ is constant555For clarity, we reuse the species labels and orientation of the gradient used in the present work. This should be kept in mind when comparing with the work of Chen et al where the gradient is along x and the species are labelled $F$ for the fuel (here, $A$), $S$ for the inert fluid (here, $F$) and $P$ for the reaction product (here, $B$)., i.e. $c_{A}(x,y)=\rho\frac{y}{L_{y}}$ and $c_{F}(x,y)=\rho\frac{L_{y}-y}{L_{y}}$ in the absence of chemical reaction, and the interaction parameters $\epsilon$ between the $C$ sphere and both solvent species is identical, $\epsilon_{C,A}=\epsilon_{C,F}=\epsilon_{C,B}$, so that there is no chemical gradient force on the $C$ sphere (explicitly, $\Lambda_{C,A}=\Lambda_{C,B}=\Lambda_{C,F}$).

The main finding of Chen et al are

The average orientation of the nanomotor is toward $y=0$, even though by a small amount (see inset of Fig. 4(c) of Ref. [13]). 2. 2.

The average trajectory of the nanomotor is toward $y=L_{y}$ (see Fig. 7(a) of Ref. [13]).

We will show here that these results depend strongly on the choice of parameters for the interaction between the colloid surfaces and the solvent via three observations: the distribution of trajectories $y(t)$, the distribution of angles $P(\theta)$, and the distribution of position $P(y)$.

We have performed stochastic simulations with $\epsilon_{\kappa,A}=\epsilon_{\kappa,F}=\epsilon_{C,\alpha}=1$ (for all solvent species $\alpha$ and colloid species $\kappa$) and $\epsilon_{N,B}=0.25$, $L_{y}=30$, $L_{z}=15$. A harmonic wall with spring constant $k_{\mathrm{wall}}=10$ prevents the exit of the colloid and is turned on for $y<2$ or $y>\L_{y}-2$. The nanomotors were started in the $x$ direction, as in Ref. [13], to avoid an a priori bias. The results are displayed in Fig. 11. The distribution of $y(t)$ trajectories in the upper panel shows a large spread, from which we cannot conclude of a dominant chemotactic behaviour. Examining $P(y)$, we see however that there is an average accumulation of particles close to $y=0$. The angular distribution $P(\theta)$ does show a bias toward $y=0$ (equivalently $\theta=\pi$), as in Ref. [13].

Until here, there was no chemically-induced force on the $C$ sphere. To assess the importance of this choice, we perform another set of simulations $\epsilon_{C,B}=\epsilon_{N,B}$ that are shown in Fig. 12. This parametrisation is the one used for our main results of sections 6.1, 6.2, and 6.3. From Eqs. (31)-(34), we know that the torque on the nanomotor will be influenced by this choice. This is reflected in the angular distribution $P(\theta)$ in Fig. 12 that now displays a strong orientation toward $\theta=0$ (i.e. the nanomotor is oriented toward $y=L_{y}$). As a result, there is no competition between a downward-facing nanomotor and the greater velocity for upward-facing nanomotor orientation that is observed in Ref. [13] and the distribution of $y(t)$ trajectories is narrow, showing a strong chemotactic behaviour. This is confirmed by the distribution of positions $P(y)$.

To conclude this comparison, we have observed competing effects of orientation and propulsion, leading to different chemotactic behaviours. The distribution of angles, for instance, depends stronly on the parameters chosen for the colloid-fluid interaction: For the interaction choice of Chen et al, $\epsilon_{C,A}=\epsilon_{C,F}=\epsilon_{C,B}$, we also find an average orientation toward $y=0$, but for the other choice that we used the average orientation is toward $y=L_{y}$. The overall chemotactic behaviour will however also depend on the diffusive and propulsive properties of the dimer.

Although further research on the experimental characterisation of the surface properties of the nanomotors is needed, our stochastic model provides an interesting tool to relate these properties to the effective chemotactic behaviour of nanomotors.

7 Conclusions

We have proposed a particle-based simulation setup, based on Multiparticle Collision Dynamics and Molecular Dynamics, to study the chemotactic motion of passive and active colloids. The concentration field that drives this motion is sustained by the flowing input of the cell, via two inlets, as is done in the experiment of Ref. [9].

We then constructed an approximate solution for the chemical concentration profile in the microfluidic channel in the presence of an active colloid and built a stochastic model with a microscopic expression for the systematic force. This model provides a sound basis for the microscopic origin of chemotaxis, i.e. it explains why the colloids move towards higher or lower values of the coordinate $y$. Upon extending this model to a two-sphere dimer nanomotor, we also gain an understanding of why the nanomotor changes its orientation via a systematic torque. We have thus improved on the formal understanding of the chemotactic motion in the microfluidic channel.

Our stochastic model, adapted to the setup with a constant concentration gradient by Chen et al reproduces the observation that the nanomotor tends to orient opposite to the gradient. While Chen et al observed a net positive chemotaxis, this is not the case in the present work. The reason is that the interplay between orientation and propulsion depends on the geometry of the motor and on the choice of parameters. This reasoning is supported by changing the surface interaction also for the $C$ bead, resulting in a cooperation of orientation and propulsion leading to positive chemotaxis. The specific change in $\epsilon_{C,B}$, that is not used in the simulation literature but introduced here in subsection 6.3, allowed us to probe a regime of well-defined positive chemotaxis and demonstrates that the class of simulation models introduced by Rückner and Kapral [27] possesses a very rich phenomenology.

It is interesting to note that our stochastic model requires only the input of the surface interaction parameters $\Lambda_{\kappa,\alpha}$ and of the diffusion coefficients of the dimer. In principle, it could be extended to other forms of the interaction potential or rely on a characterisation of $\Lambda_{\kappa,\alpha}$ that does not reveal the full functional form of this potential, and also to other motor geometries. While we have chosen $\Lambda_{N,\alpha}=\Lambda_{C,\alpha}$ for all solvent species $\alpha$ for this first investigation of the model, there is room for possible qualitative changes with other choices, as we have witnessed for the constant gradient configuration.

The overall speed of the motor, as it was the case in earlier simulation studies for dimer nanomotors, only depends on $\Lambda_{N,A}-\Lambda_{N,B}$ [27]. The torque, on the other hand, depends also on $\Lambda_{C,A}-\Lambda_{C,B}$. A difference in these quantities, and thus in the surface properties of both sides of the nanomotor, could be revealed in a controlled manner in experiments and this reinforces the importance of the microfluidic channel configuration for chemotactic studies.

The present work can be extended to other types of motors, notably the Janus nanomotors that are used in Ref. [9] and for which colloidal assemblies have already been used in mesoscopic simulations [41]. This line of research is promising to test in silico the behaviour of different motor geometries. Using models suitable for the chemo-mechanics of enzymes, at a mesoscopic level [42, 43, 44], could provide very fruitful advances for understanding the recent works on enzyme chemotaxis [11, 12], especially given the fact that multiple inlets microfluidic devices originate from studies on bacterial chemotaxis [45] and have been used for the enzyme studies [12, 11].

Acknowledgements

The authors thank Raymond Kapral for interesting discussions and Samuel Sanchez for communicating the height of the cell in Ref. [9]. PdB is a postdoctoral fellow of the Research Foundation-Flanders (FWO).

Appendix A Computational reproducibility

In this appendix, we review how the present work can be reproduced. The software and parameter files are all available publicly under open-source licences. We have prepared supplementary material that contains the relevant parameter files for the mesoscopic simulations and the code for data analysis and archived them with Zenodo666https://zenodo.org/, available as Ref. [46].

All mesoscopic simulations are performed using the open-source software package RMPCDMD [29, 30] for the simulations of passive and active colloids, developed by the authors with Mu-Jie Huang and Peter Colberg. The output of RMPCDMD consists of H5MD [47] files that contain the full trajectory for the colloids, the thermodynamic observables and the correlation functions (velocity autocorrelation functions and mean-squared displacement).

All stochastic simulations are performed using Python, NumPy and Cython in a Jupyter777http://jupyter.org/ notebook. The analysis of both types of simulations and the execution of the stochastic model simulations are done in the notebook colloidal_chemotaxis.ipynb, except for the constant gradient model that is implemented in a separate Cython module stochastic_nanomotor.pyx and driver program run_cg_nm.py. The equilibrium simulations for the dimer are analysed in the notebook diffusion.ipynb.

The references for the software are the following: NumPy [48] is used for all numerical work in the analysis of the mesoscopic model and overall for the stochastic model, SciPy [49] is used for computing the erf function and for numerical integration, matplotlib [50] to generate the figures, Mayavi [51] for Fig. 2, h5py [52] to read simulation data, Cython [53] to accelerate the nanomotor stochastic simulations, gfortran [54] to build the RMPCDMD code.