Drag force on a microrobot propelled through blood

Problem statement and basic equations

Consider a rigid spherical robot driven by an external magnetic force along the –z direction in an infinite suspension of RBCs modeled as N RBCs suspended in a triply periodic suspension. The magnetic thrust force F_m imparted by a magnetic field B existing in z direction only, written as:

where a_s is the radius of the robot and M is the magnetization. We assume that the density of the robot is the same as that of the ambient fluid (ρ), and do not consider sedimentation effects. The ambient fluid is a Newtonian liquid of viscosity μ. Given the small size of the robot, the surrounding flow is Stokes flow. As the effect of inertia is negligible, the Stokes’ drag force and the magnetic thrust F_m are balanced. The speed of the robot (v_N) in a Newtonian fluid of viscosity μ can be derived as:

Note that this formula cannot be used when RBCs are present because it is assumed that the ambient fluid is a homogeneous fluid of viscosity μ. Each RBC is modeled as a capsule with a hyper-elastic membrane. We assume that the internal and external RBC fluids have the same density ρ and viscosity μ. RBC deformation becomes more evident as the magnetic thrust of the robot increases. The ratio between the magnetic force and elastic force of the RBC membrane, Γ, is an important dimensionless magnetoelastic number that governs RBC deformation and is defined as:

where G_s is the shear elastic modulus of RBC membrane and a_RBC is the characteristic radius of an RBC defined as \(\root 3 \of {\frac{3{V}_{{{{{{{{\rm{R}}}}}}}}}}{4\pi }}\), where V_R is the volume of an RBC. Another important dimensionless number is ϵ, which is defined as:

which determines the size effect of the robot relative to that of RBCs.

When a robot navigates in an infinite RBC suspension (RBCs are assumed to be neutrally buoyant), i.e., a suspension with N RBCs in a triply periodic flow field, the the velocity at a given point x is given by the boundary integral formulation⁶⁹:

where A_R is the RBC surface and A_s is the robot surface. q_R and q_s are the force densities on the RBC membrane and surface of the robot, respectively. J is the Green’s function, given by:

where δ_ij is the Kronecker delta, r = x − y is the vector from the source point y to the observation point x, and r = ∣r∣. J ^E is Green’s function for the triply periodic lattice, which is solved using the Ewald summation technique, written as⁷⁰:

where l_x, l_y, and l_z are the lattice sizes in the x, y and z directions, respectively, and \(\xi = {\pi \,}^{\frac{1}{2}}{({l}_{x}{l}_{y}{l}_{z})}^{-\frac{1}{3}}\). γ is the index of periodic boxes and κ is the index of reciprocal vectors. k is the reciprocal vector and k = ∣k∣. n_x, n_y, n_z are the integer numbers of the indices of the periodic boxes. To avoid RBC-RBC and RBC-robot overlaps, a very short-range repulsive force F_rep is employed, which acts on the observation Gauss point \({{{{{{{{\bf{y}}}}}}}}}_{j}^{{{{{{{{\rm{G}}}}}}}}}\) and is expressed as:

where k_c is a spring constant and d_G is defined as:

\({{{{{{{{\bf{x}}}}}}}}}_{i}^{{{{{{{{\rm{G}}}}}}}}}\) are the nearby Gauss points of distinct RBCs. In this study, k_c is set to 0.04G_s and the repulsive force is considered only when the distance between two Gauss points on different RBCs is smaller than 0.18a_RBC.

The RBC membrane is assumed to be an isotropic and hyperelastic material that follows the Skalak law⁴⁵. The surface deformation gradient tensor D_s is given by:

where X and x are the membrane material points of the reference and deformed states, respectively. Local deformation of the RBC membrane is described by the Green-Lagrange strain tensor:

where I_s is the tangential projection operator. Two invariants of the in-plane strain tensor E are:

where λ₁ and λ₂ are the principal stretch ratios. The Cauchy stress tensor T is:

where S = λ₁λ₂ is the area dilation ratio. Invoking the two-dimensional constitutive law of Skalak, the elastic strain energy per unit area w_s can be expressed as:

where C is the area dilation modulus (10 in this paper)⁵².

Numerical method

The numerical method developed by Walter et al.⁴⁵ was applied to couple the BEM of fluid mechanics to the FEM of solid mechanics. We used this method in previous papers^48,49,50,51. A single RBC membrane is discretized using an unstructured triangular mesh produced via recursive subdivision of an icosahedron. In solid mechanics, the weak form of the equilibrium condition of the RBC is:

which is solved by the FEM, where \(\hat{{{{{{{{\bf{u}}}}}}}}}\) and \(\hat{{{{{{{{\boldsymbol{\epsilon }}}}}}}}}\) are the virtual displacement and the virtual strain, respectively. In fluid mechanics, the BEM combined with multipole expansion is applied to solve Equation (11). The surface integrals of the BEM are computed using the Gaussian quadrature approach. Polar coordinates are utilized to eliminate the 1/r singularity^53,71. The velocity of the RBC membrane at point x is given by the kinematic condition:

The position is updated using the explicit second-order Runge-Kutta method. The implementation of a volume constraint prevents even a modest volume error⁷¹.

When resolving the velocity disruptions caused by near-field forces, a fine mesh (small element size) is vital; however, the velocity caused by far-field forces can be coarse-grained when using the adaptive mesh method. In the simulations, when a certain RBC (source set) is distant from another RBC or the robot (including the observation point), a patch of 16 elements is used to replace the triangular elements. Thus, an RBC is constructed from 80 large patches. In contrast, an RBC has 1,280 triangular elements when a fine mesh is required. The force density on a patch with N_e fine elements (q_c) is given by:

where \({q}_{{{{{{{{\rm{f}}}}}}}}}^{i}\) is the force density exerted on fine mesh.

To enhance the stability of simulations, the second-layer computation of the Ewald summation was introduced. Unfortunately, this greatly increased the computational burden. To reduce this load while retaining accuracy, the boundary integral of q_R in Eq. (11) is expanded in terms of moments^72,73:

where V = l_xl_yl_z, and H. O. T. are the higher order terms. The force and torque terms disappear due to the force- and torque-free conditions. The stresslet H^m is as follows:

where \(\hat{{{{{{{{\bf{r}}}}}}}}}\) is the vector from the center of a certain RBC to a defined point on the membrane of the same RBC. The propagators are:

We used the stresslet as a substitute for an RBC during second-layer Ewald summation; this significantly accelerated the computation. The nomenclature and parameter ranges used in this study are summarized in Supplementary Table 3.

Experimental method and materials

To validate the simulation, we conducted experiments involving the sedimentation of spherical particles in an RBC suspension. Similar to the magnetic force used in the simulation, gravity force is also considered a body force, and it can serve as a driving force for the particles. In this context, we denote Γ as \(\it {{\varGamma }}_{\exp }\) and define it as follows:

Here, ρ_p and g are the particle density and the gravitational acceleration. \(\it {{\varGamma }}_{\exp }\) still represents the ratio of driving force and elastic force acting on the RBC membrane, although the driving force has shifted from a magnetic force to a gravity force.

The preparation of the RBC suspension, following the same procedure as described in our previous research^74,75,76, is briefly outlined here. Pig blood was bought from Tokyo Shibaura Organ Corporation. The RBCs were separated from the bulk pig blood through centrifugation at 1500 rpm for 20 minutes (model 3220 Micro Hematocrit Centrifuge, Kubota). Subsequently, plasma and buffy coat were removed by aspiration. The washing and centrifuging with physiological saline (PS) were repeated twice. This process of washing and centrifugation with physiological saline (PS) was repeated twice. The washed RBCs were then diluted with dextran 40 (DX-40, Otsuka Medicine)

The observation chamber was assembled using two cover glasses (NEO Micro cover glass; Thickness No. 1; 24 × 60 mm, Matsunami, Tokyo, Japan) and five layers of double-sided tape. The height of the chamber was approximately 1500 μm. Initially, the top cover glass, the observation chamber without the top cover glass, and iron particles (d50 = 5.0 μm, type REZL, Jiangsu Tianyi Ultrafine Metal Powder Co. Ltd., China) were immersed in bovine serum albumin (BSA, Stock Solution 130091-376, Miltenyi Biotech., Germany) for an hour to prevent particle aggregation and adhesion to pig RBCs or the inner surface of the observation chamber^77,78. Subsequently, a mixture of the particle suspension and the RBC suspension was prepared to create the test suspension containing 25% RBCs (25% hematocrit), 15% RBCs, and 5% RBCs. The test suspension was injected into the observation chamber after BSA was aspirated. Finally, the observation chamber was sealed via the top cover glass. For the control group, the particle suspension and pure dextran 40 were mixed in the same proportion as in the test suspension. Consequently, all components in the test suspension and control suspension, except for RBCs, are identical.

Sample examination was performed using an inverted microscope (Olympus IX71; Olympus, Tokyo, Japan) equipped with a 40× magnification objective, as outlined in Fig. 5b. Initially, the test suspension or control suspension was placed on the microscope stage for 30 min, allowing the complete particle settlement (Fig. 5c ① ②). The microscope’s focus was then finely adjusted to yield a sharp image of the particles. Subsequently, the observation chamber was inverted and left in that position for an additional 30 minutes to direct particle settlement onto the coverslip (Fig. 5c ③). Finally, the observation chamber was inverted again after video recordings were captured using a high-speed camera operating at a frame rate of 50 fps (Fig. 5c ④ ⑤). For our analysis, only particles with diameters ranging from 4 μm to 6 μm were considered. A snapshot displaying the sedimentation at the chamber’s bottom is shown in Fig. 5d. Sedimentation time could be calculated from the first frame that exhibits the stationary observation chamber to the frame where a distinct spherical iron particle becomes visible. The sedimentation velocity of particles can be calculated by dividing the depth of the chamber by the sedimentation time of the particles. Since the driving force, gravity, is constant, this sedimentation velocity can be compared to that without RBCs to obtain the relative resistance coefficient of a microrobot (\({C}_{{{{{{{{\rm{r}}}}}}}}}^{* }\))

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.