Modeling of 3D Cilium Mechanics in COMSOL® Using Beam Physics with Linear Extrusion Coupling

Y. Shen1, L. G. Woodhams1, P. V. Bayly1
1Washington University in St. Louis, St. Louis, Missouri, USA
Published in 2020

Cilia are slender micro-organelles (200 nm diameter, 10 µm long) that generate propagating waves to propel cells or move fluid. The cytoskeletal structure of the cilium (the axoneme) consists of 9 outer microtubule doublets and 2 central microtubule singlets. Outer doublets are connected by active motor proteins (dyneins) and passive components (nexin links). The structure of the cilium is well known by electron microscopy, but the mechanism by which cilia oscillate remains uncertain. By building a 3D finite-element model of the axoneme, we can investigate the theory behind cilia motion and evaluate the effects of system parameters on stability. In fact, using this model, we can estimate coupling parameters that we are not able to measure experimentally in these nanoscale structures.

While the actual structure of a cilium is more intricate, we can model the biological system with beam structures, including a four-outer-doublet system (4-beam model) as a theoretical baseline and a six-outer-doublet with a central-singlet system (6+1 beam model) as a more realistic model. COMSOL Multiphysics® beam system module is used to build the model. Both of the models are further simplified by limiting loading to two active pairs connected by passive links. Linear extrusion coupling operators are used to calculate average tangential vectors for the active pairs and map relative positions of different beams. For the active pairs, “follower” loads modeling dynein forces are applied along a direction that is the average of the tangent vectors of the beams within the active pairs. For the passive pairs, the inter-doublet forces are calculated based on nonlinear viscoelastic spring equations with displacements probed using the extrusion coupling operator. Reaction moments for the dynein arms are applied by the cross product of relative displacement and “follower” load. The model is solved using the time dependent solver for a characteristic time period. Presence of oscillatory behavior is determined through 1D point plots of tip displacement, and characteristic beating shape is observed through 3D plots and animations.

Three predominant behaviors were observed through parametric sweeps: (i) stable equilibria; (ii) flutter; and (iii) divergence. In the stable regime, dynein forces are not enough to drive the system to motion, leaving the system in a static equilibrium with no bending. In the flutter region, the system oscillates periodically; this behavior corresponds to complex conjugate pairs of eigenvalues with positive real parts in the least stable mode. In the divergence region, bending deformations grow without oscillation; the least stable mode has an eigenvalue with a positive real part. Using these simplified 3D cilium models, we observed how varying parameters can cause transitions between these three behaviors. The model also opens several questions for future work, most notably: Does one of these instabilities (flutter or divergence) or a combination of both underlie the true mechanism behind the wavelike oscillation of cilia?