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

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 nanomotors 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.

---

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,¹ polystyrene² and silica³ Janus particles, dimers⁴). 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–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.⁸ 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.¹⁰ 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.¹¹ 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.¹³ 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.⁹ 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 build 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.¹³ 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 mass m_i, 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 τ 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,¹⁶ or by Nikoubashman et al. to study the flow of colloids in the presence of obstacles,¹⁷ 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

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 MPCDwhere N_MD is the number of MD steps between successive MPCD collisions.

At fixed time intervals τ, the fluid particles are sorted in a lattice of cubic cells of side a, whose origin is shifted randomly to ensure Galilean invariance,²⁰ and their velocities are collided cell-wise according to

v_i′ = v_ξ + Ω_ξ(v_i − v_ξ) (4)

where the prime denotes the post-collision value, ξ is the cell containing particle i, Ω_ξ a rotation operator of angle Θ in [Doublestruck R]³ around a randomly chosen axis and v_ξ is the centre-of-mass velocity of the cell. The transport properties of a MPCD fluid can be computed analytically, see ref. 21, 22 and references therein. The simulation parameters are given in 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	ρ	10
Cell size	a	1
MPCD time step	τ	0.5
MD steps per MPCD step	N MD	50
Cell dimensions	L x , Ly, Lz	90, 60, 15
Buffer length	L buffer	20
Temperature	k B T	1/3
Acceleration	g	0.001
MPCD collision angle	Θ	2.6
Reaction probability	p	1
Colloid radius	σ	3
Unit interaction energy	ε A	1

---

One technique to obtain a Poiseuille flow is to impose a constant acceleration field in the simulation cell,²³ 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 ref. 19 and 24, 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.²⁵ 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 acceleration field g. In the case of dimer nanomotors, the distance between the two spheres is held constant using the RATTLE algorithm.²⁶

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

where ε_ij denotes the strength of the potential and σ_ij denotes 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 ε_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. ε_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 ε_A and time . As the interaction between all colloids and solvent particles of species A is equal, the first subscript to ε_κ,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

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²⁷ or in the bulk.²⁸ At the surface of a catalytic colloid C, the reaction

A + C → B + C (7)

converts fluid particles of type A (the fuel) to fluid particles of type B (the product). The reaction occurs when the fluid particles cross the interaction region of the colloid and is executed with a probability p ∈ [0,1], when the fluid particles exit in the interaction region to avoid any discontinuity in the energy.²⁷

All simulations were performed using 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.⁹ 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. 1. 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 or 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_buffer. In the following text of the paper, the trajectories from the mesoscopic simulations are shifted in x by −L_buffer to match the coordinate system of the stochastic model.

Fig. 1 Schematic description of the simulation cell for chemotaxis in a flow. There are two inlets on the left. In the buffer (leftmost region), solvent particles are set to species A for 0 < y < L_y/2 and F else. Lennard-Jones 9-3 potentials confine the colloids close to z = L_z/2. The colloid are initially placed in the F inlet and constrained to a fixed y and z track. This constraint is released when x > σ.

Fig. 2 Simulation snapshot for a dimer nanomotor simulation. The confining plates are made of transparent blue. The dimer is made of a catalytic (red) bead and of a non-catalytic (blue) bead and its centre-of-mass trajectory is shown as a white line. Here, the nanomotor is attracted toward the front of the cell (lower values of y) and displays a corresponding reorientation. The instantaneous concentration field of A solvent particles around z = L_z/2 is shown in pseudocolor.

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_flow and average velocity v_av = 2/3 v_flow. We report the characteristic numbers of the flow in Table 2. For the work of Baraban et al.,⁹ we use the values found in the article for the flow rate (140 μL per hour), the width (200 μm per inlet, there are three inlets) and temperature (300 K). The height of the channel was confirmed by email to be 30 μm. For the properties of water, we use reference data from the NIST³¹ (retrieved January 11, 2017: density ρ_water = 996.56 kg m⁻³, and viscosity η_water = 8.54 × 10⁻⁴ Pa s) and take the diffusion coefficient D_exp = 2 × 10⁻⁹ m² s⁻¹ that is a reasonable approximation for both water³² and hydrogen peroxide.³³ For the MPCD fluid properties, we refer to the review by Kapral.²² The viscosity is

and the self-diffusion coefficient isWe have verified that the velocity profile obeys the theoretical prediction²⁴

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 (Ly)	60	600 μm
Height (Lz)	15	30 μm
Average flow velocity vav	0.063	2.16 mm s^−1
Maximum flow velocity vmax	0.095	3.24 mm s^−1
Pe	14	32
Re	0.48	0.076

---

We have aimed for a similar fluid regime with respect to the experiment, that is a high Péclet number (Pe) flow, as our estimate for the concentration profile relies on Pe ≫ 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 Re should remain moderate for the laminar regime to hold.

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

Pe = L_zv_av/D,

Re = v_flowL_z/η,

and
where is the speed of sound.¹⁶ The maximum Mach number is obtained for the maximum velocity of the flow and is Ma ≈ 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 and . This restriction is lifted when the x position of the colloid (centre of mass) has passed L_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_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_α(r) for the different chemical species in the cell and the resulting chemotactic force on the colloids.

We approximate the evolution of c_α(r) by a 1D diffusion equation where the spatial direction of the flow is proportional to time, i.e. x = v_flowt, and the diffusion process acts in 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_α(r) in the absence of catalytic particles.

Fixing z = L_z/2 and identifying the time and the x coordinate results in

The separate inlets for the A and F chemical species translate into the initial condition c_A(0,y) = c₀Θ(L_y/2 − y), c_B(0,y) = 0 and c_F(0,y) = c₀Θ(y − L_y/2) for x = 0 (equivalently t = 0). The solution to eqn (11) is

c_F(x,y) = c₀ − c_A(x,y) (13)

c_B(x,y) = 0 (14)

where 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,in the low Mach number conditions here.

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

and are represented in Fig. 3.

Fig. 3 The spherical coordinates centered on the colloid.

The presence of a catalytic colloid is taken into account in eqn (11) 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σ, 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³⁵ and rederived later by other authors.^36,37 The RBC at radius R in the absence of external gradient, for the species A, is

Rk_D∂_rc_A = k₀c_A, (17)

where k_D = 4πRD is the diffusion-limited rate constant and is the chemical rate.³⁸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,²⁷ we sum the interaction potential V_κ,α between a colloid of species κ and the fluid particles of species α over the interfacial region to obtain the potential energy. Differentiation with respect to the colloid's coordinate gives the force on the colloid: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_α by the Boltzmann weight.²⁷Eqn (18) then reduces towhere β = (k_BT)⁻¹ and 1⃑_r is the unit vector pointing in the radial direction. Integrating by parts yieldswhere we have definedWithin this framework, it is sufficient to compute surface integrals over the colloids to compute the force acting on them. The use of eqn (20) 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 (12)–(14) together with eqn (20).Locally around the colloid, we use a first-order expansion of c_α in the lateral direction.where λ = ∂_yc_A(x,y). The explicit dependence of F_N,y and λ on the coordinates is omitted for brevity.

4.3 Single active colloid

To take into account the catalytic activity of a colloid, the diffusion eqn (11) is solved with the radiation boundary condition (RBC). Eqn (17) expresses the flux that originates from the chemical activity of the colloid, which 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 eqn (17) and of the diffusive flux†.

Rk_D∂_rc_A = k₀c_A + Rk_Dλ cos θ (25)

We make the following ansatz for the concentration field c_A:Eqn (26) is a solution of the diffusion equation and matches the boundary condition (25). 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 of the external gradient. We obtain the coefficients by inserting eqn (26) in eqn (25):where we have used for c₀ the solution (12). The value of λ is obtained by differentiating eqn (12) with respect to y.

It is important to mention that (i) in the absence of a gradient, the coefficients in eqn (27) lead to the standard solution of a spherical sink and (ii) in the absence of chemical activity, the solution (12) 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 that 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

To compute the force in eqn (20), we set Λ_C,F = Λ_C,A. Using eqn (15), we find that the chemotactic force on the active colloid C is

4.4 Dimer nanomotor

A dimer nanomotor is made of two spheres linked rigidly, with 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 on 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 eqn (27) and (29). 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

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

Fig. 4 Coordinate system for the dimer nanomotor.

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

The torque on the dimer is

The motion of the dimer is studied via its centre-of-mass position in the x–y plane and its inclination ϕ 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.³⁹where v_flow is the average velocity due to the flow at z = L_z/2, ξ_x and ξ_y are normal distributed white noise, D is the diffusion coefficient of the colloid and γ its friction. For the single spheres, we use γ = 4πησ, the slip hydrodynamics friction, and D = k_BT/γ.

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

where the projected forces F_∥ and F_⊥ areand [scr T, script letter T] is the torque defined in eqn (32). The rotation operations in eqn (36) and (37) 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 L⃑ = (32, 32, 15) and with no chemical reaction, all other simulation parameters are taken equal. The cartesian and angular velocity distributions for the dimer are shown in Fig. 5.

Fig. 5 The cartesian and angular velocity distributions for the dimer in the equilibrium simulations. The blue bars form a histogram over the simulation data and the full orange lines are the Boltzmann distributions at k_BT = 1/3.

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

where ζ is ∥ or ⊥, following ref. 39. The projected velocities are defined bywith û the unit vector in the direction joining the two spheres. The angle ϕ measures the deviation of û in the x–y plane with respect to the unit vector 1⃑_x. The results of these simulations are D_∥ =2.0 × 10⁻³ and D_⊥ = 1.5 × 10⁻³. The rotational motion is characterised in the same fashion, using the velocity autocorrelation of angle ϕ: 〈[small phi, Greek, dot above](τ)[small phi, Greek, dot above](0)〉 with the result D_r = 1.4 × 10⁻⁴. We determine the friction coefficient using the fluctuation–dissipation relation The rotational time D_r⁻¹ ≈ 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 ε_A ≠ ε_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, ε_N,F = ε_A so there is no chemotactic behaviour.

Fig. 6 Example mesoscopic simulation of a passive sphere with ε_N,F = 1. The black line denotes the trajectory of the colloid that starts on the left at the entrance of the cell and follows the flow to the right. The motion of the colloid is constrained to y = L_y/2 + y_shift until the release point x = σ (denoted by a black circle), after which it is influenced by the chemotactic force. The pseudocolor background indicates the strength of the gradient λ(x,y) = ∂_yc_A(x,y).

Upon changing ε_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 the passive diffusiophoresis.

The results for passive spheres are summarised in Fig. 7 for all chosen values of ε_N,F. For the central panel ε_N,F = ε_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 ε_N,F < ε_A, we have Λ_N,A − Λ_N,F < 0 and we expect from eqn (24) a lateral force whose sign is opposite to gradient λ. As c_A(x,y) decreases for increasing y, λ < 0 and the chemotactic force is upward, as confirmed by the panels ε_N,F = 0.25 and ε_N,F = 0.50. Conversely, for ε_N,F > ε_A, the chemotactic drift is negative. This is shown for ε_N,F = 2 and ε_N,F = 4 in Fig. 7.

Fig. 7 Ensemble trajectories for the passive sphere simulations. The orange (blue) line is for the average position of the colloid in the mesoscopic (stochastic) simulations and the corresponding filled area indicates ± one standard deviation across realisations. There are 16 realisations for every set of parameters and simulation type.

Fig. 7 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 Λ_N,A − Λ_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 that for the passive sphere but now ε_C,F = ε_A and the mechanism that creates a systematic lateral motion for the passive sphere, that is proportional to Λ_C,A − Λ_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 Λ_C,B − Λ_C,A, creates an original combination of passive and active diffusiophoresis. Using eqn (29) and observing that c₂ > 0 (as λ < 0 everywhere in the simulation), we expect a downward chemotaxis for Λ_C,B − Λ_C,A < 0 (i.e. ε_C,B < ε_A) and an upward chemotaxis for Λ_C,B − Λ_C,A < 0 (i.e. ε_C,B > ε_A). We explore the effect of changing ε_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.

Fig. 8 Ensemble trajectories for the active sphere simulations. Organisation as in Fig. 7.

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 r⃑_com and the orientation ϕ. The release from the injection track occurs here as soon as both spheres have exited the buffer region and the initial orientation is ϕ = 0.

We have chosen, for the simulation parameters, that ε_C,α = ε_N,α, for all solvent species α. This is at variance with simulations of the prototypical dimer nanomotor^4,13,27,40 where ε_C,A = ε_C,B = ε_A ≠ ε_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 ϕ 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 ϕ are reported in Fig. 9 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 eqn (35) and (36) as both the directions in y and in ϕ 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 eqn (27), (29) and (31). 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.

Fig. 9 Ensemble trajectories for the dimer nanomotor simulations. The centre-of-mass position is used. Organisation as in Fig. 7.

Fig. 10 Ensemble trajectories of the dimer nanomotor simulations. Here, the orientation ϕ with respect to the x axis is shown. Organisation as in Fig. 7.

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. designed a cell in which the gradient of the chemical species A and F is constant,‡i.e. and in the absence of chemical reaction, and the interaction parameters ε between the C sphere and both solvent species are identical, ε_C,A = ε_C,F = ε_C,B, so that there is no chemical gradient force on the C sphere (explicitly, Λ_C,A = Λ_C,B = Λ_C,F).

The main findings of Chen et al. are

• The average orientation of the nanomotor is toward y = 0, even though by a small amount (see the inset of Fig. 4(c) of ref. 13).

• 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(θ), and the distribution of position P(y).

We have performed stochastic simulations with ε_κ,A = ε_κ,F = ε_C,α = 1 (for all solvent species α and colloid species κ) and ε_N,B = 0.25, L_y = 30, L_z = 15. A harmonic wall with spring constant k_wall = 10 prevents the exit of the colloid and is turned on for y < 2 or y > Ł_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 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(θ) does show a bias toward y = 0 (equivalently θ = π), as in ref. 13.

Fig. 11 Simulations of the constant-gradient setup with ε_N,B = 0.25. Averages are performed over 40 realisations for 2.5 × 10⁴ < t < 10 × 10⁴ (neglecting the start of the simulations for about 3 typical rotation times 1/D_r). From top to bottom, the panels show the area between 〈y(t)〉 − σ_y(t) and 〈y(t)〉 + σ_y(t), the distribution of angles θ with respect to the y axis and the distribution of y.

Until now, there was no chemically-induced force on the C sphere. To assess the importance of this choice, we perform another set of simulations ε_C,B = ε_N,B that are shown in Fig. 12. This parametrisation is the one used for our main results of Sections 6.1–6.3. From eqn (29)–(32), we know that the torque on the nanomotor will be influenced by this choice. This is reflected in the angular distribution P(θ) in Fig. 12 that now displays a strong orientation toward θ = 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).

Fig. 12 Simulations of the constant-gradient setup with ε_N,B = 0.25 and ε_C,B = 0.25. Averages are performed over 20 realisations for 2.5 × 10⁴ < t < 10 × 10⁴. Panels as in Fig. 11.

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 strongly on the parameters chosen for the colloid–fluid interaction: for the interaction choice of Chen et al., ε_C,A = ε_C,F = ε_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 described in 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 ε_C,B, that is not used in the simulation literature but introduced here in Section 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²⁷ possesses a very rich phenomenology.

It is interesting to note that our stochastic model requires only the input of the surface interaction parameters Λ_κ,α 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 Λ_κ,α that does not reveal the full functional form of this potential, and also to other motor geometries. While we have chosen Λ_N,α = Λ_C,α for all solvent species α 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 Λ_N,A − Λ_N,B.²⁷ The torque, on the other hand, depends also on Λ_C,A − Λ_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.⁴¹ 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–44 could provide very fruitful advances for understanding the recent works on enzyme chemotaxis,^11,12 especially given the fact that multiple inlet microfluidic devices originate from studies on bacterial chemotaxis⁴⁵ and have been used for the enzyme studies.^11,12

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 Zenodo,§ 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⁴⁷ 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 Jupyter¶ 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⁴⁸ is used for all numerical work in the analysis of the mesoscopic model and overall for the stochastic model, SciPy⁴⁹ is used for computing the erf function and for numerical integration, matplotlib⁵⁰ to generate the figures, Mayavi⁵¹ for Fig. 2, h5py⁵² to read simulation data, Cython⁵³ to accelerate the nanomotor stochastic simulations, gfortran⁵⁴ to build the RMPCDMD code.

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).

References

1. W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. S. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert and V. H. Crespi, J. Am. Chem. Soc., 2004, 126, 13424–13431.
2. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 48102.
3. H. Ke, S. Ye, R. L. Carroll and K. Showalter, J. Phys. Chem. A, 2010, 114, 5462–5467.
4. L. F. Valadares, Y.-G. Tao, N. S. Zacharia, V. Kitaev, F. Galembeck, R. Kapral and G. A. Ozin, Small, 2010, 6, 565–572.
5. R. Kapral, J. Chem. Phys., 2013, 138, 020901.
6. J. Wang, Nanomachines, Wiley-VCH, Weinheim, Germany, 2013.
7. S. Ebbens, Curr. Opin. Colloid Interface Sci., 2016, 21, 14–23.
8. Y. Hong, N. M. K. Blackman, N. D. Kopp, A. Sen and D. Velegol, Phys. Rev. Lett., 2007, 99, 178103.
9. L. Baraban, S. M. Harazim, S. Sanchez and O. G. Schmidt, Angew. Chem., Int. Ed., 2013, 52, 5552–5556.
10. H. C. Berg, Annu. Rev. Biophys. Bioeng., 1975, 4, 119–136.
11. S. Sengupta, K. K. Dey, H. S. Muddana, T. Tabouillot, M. E. Ibele, P. J. Butler and A. Sen, J. Am. Chem. Soc., 2013, 135, 1406–1414.
12. K. K. Dey, S. Das, M. F. Poyton, S. Sengupta, P. J. Butler, P. S. Cremer and A. Sen, ACS Nano, 2014, 8, 11941–11949.
13. J.-X. Chen, Y.-G. Chen and Y.-Q. Ma, Soft Matter, 2016, 12, 1876–1883.
14. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605.
15. A. Malevanets and R. Kapral, J. Chem. Phys., 2000, 112, 7260.
16. C. Prohm, M. Gierlak and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 80.
17. A. Nikoubashman, C. N. Likos and G. Kahl, Soft Matter, 2013, 9, 2603–2613.
18. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987.
19. J. K. Whitmer and E. Luijten, J. Phys.: Condens. Matter, 2010, 22, 104106.
20. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 20201.
21. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, Adv. Polym. Sci., 2008, 221, 1–87.
22. R. Kapral, Adv. Chem. Phys., 2008, 140, 89.
23. J. Koplik, J. R. Banavar and J. F. Willemsen, Phys. Rev. Lett., 1988, 60, 1282–1285.
24. E. Allahyarov and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 36702.
25. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, EPL, 2001, 56, 319.
26. H. C. Andersen, J. Comput. Phys., 1983, 52, 24.
27. G. Rückner and R. Kapral, Phys. Rev. Lett., 2007, 98, 150603.
28. K. Rohlf, S. Fraser and R. Kapral, Comput. Phys. Commun., 2008, 179, 132.
29. P. de Buyl, M.-J. Huang and L. Deprez, J. Open Res. Softw., 2017, 5, 3.
30. P. de Buyl, P. H. Colberg, L. Deprez and M.-J. Huang, RMPCDMD, 2016, https://doi.org/10.5281/zenodo.198599, Software release 1.0.
31. M. M. E. W. Lemmon and D. Friend, NIST Chemistry WebBook, National Institute of Standards and Technology, Gaithersburg MD, USA, 2016.
32. J. H. Wang, J. Phys. Chem., 1965, 69, 4412.
33. D. M. H. Kern, J. Am. Chem. Soc., 1954, 76, 4208–4214.
34. R. F. Ismagilov, A. D. Stroock, P. J. A. Kenis, G. Whitesides and H. A. Stone, Appl. Phys. Lett., 2000, 76, 2376–2378.
35. F. C. Collins and G. E. Kimball, J. Coll. Sci., Imp. Univ. Tokyo, 1949, 4, 425–437.
36. M. Schell and R. Kapral, J. Chem. Phys., 1981, 75, 915–920.
37. D. T. Gillespie, J. Chem. Phys., 2009, 131, 164109.
38. K. Tucci and R. Kapral, J. Chem. Phys., 2004, 120, 8262.
39. J. K. Dhont, An Introduction to Dynamics of Colloids, Elsevier, Amsterdam, The Netherlands, 1996, vol. 2.
40. Y.-G. Tao and R. Kapral, J. Chem. Phys., 2008, 128, 4518.
41. P. de Buyl and R. Kapral, Nanoscale, 2013, 5, 1337–1344.
42. A. Cressman, Y. Togashi, A. S. Mikhailov and R. Kapral, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 050901(R).
43. C. Echeverria, Y. Togashi, A. S. Mikhailov and R. Kapral, Phys. Chem. Chem. Phys., 2011, 13, 10527–10537.
44. Y. Togashi and A. S. Mikhailov, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 8697–8702.
45. H. Mao, P. S. Cremer and M. D. Manson, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 5449–5454.
46. P. de Buyl and L. Deprez, Colloidal chemotaxis companion, 2017, https://doi.org/10.5281/zenodo.344036.
47. P. de Buyl, P. H. Colberg and F. Höfling, Comput. Phys. Commun., 2014, 185, 1546–1553.
48. S. van der Walt, S. Colbert and G. Varoquaux, Comput. Sci. Eng., 2011, 13, 22–30.
49. E. Jones, T. Oliphant and P. Peterson, et al., SciPy: Open source scientific tools for Python, 2001, http://www.scipy.org/.
50. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
51. P. Ramachandran and G. Varoquaux, Comput. Sci. Eng., 2011, 13, 40–51.
52. A. Collette, Python and HDF5, O'Reilly, 2013.
53. S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D. Seljebotn and K. Smith, Comput. Sci. Eng., 2011, 13, 31–39.
54. The GNU Fortran Compiler (4.9.2), http://https://gcc.gnu.org/onlinedocs/gcc-4.9.2/gfortran/index.html.

---

Footnotes

† The flux is multiplied by the surface of the sphere 4πR² as this is how the RBC is expressed in the literature.

‡ For 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).

§ http://https://zenodo.org/

¶ http://jupyter.org/

---