Dynamic magnetic properties of magnetosomes

Magnetosomes are found in magnetotactic bacteria, and consist of linear chains of tens of single-domain magnetic nanoparticles of Fe3O4 (magnetite) or Fe3S4 (greigite), enclosed in a lipid bilayer membrane. The particles are typically in the size range 30– 120 nm , and are practically monodisperse. Harnessing the power of biomineralisation could lead to efficient strategies for synthesising semiflexible dipolar filaments, and the development of optimum materials for applications, including in biomedicine. Brownian dynamics simulations of noninteracting magnetosomes, containing 1⩽N⩽64 ferromagnetic nanoparticles, have been used to determine static properties, and the dynamical response to a weak AC magnetic field. Results are presented for the radius of gyration Rg , the static magnetic susceptibility χ(0) , the dynamic magnetic susceptibility χ(ω) , and the effective Brownian rotation time τrot , all as functions of N. The results are compared to theoretical predictions of the flexibility, susceptibility, and rotational dynamics of such magnetic filaments.

(Some figures may appear in colour only in the online journal)

Introduction
Magnetosomes are biological structures comprised of chains of ferromagnetic or ferrimagnetic nanoparticles encapsulated in a lipid membrane [1][2][3]. These organelles are found in magnetotactic bacteria, which are thought to use the alignment of the magnetised chain with the Earth's magnetic field to aid navigation towards nutrients. The particles are iron-rich, and are most commonly formed from magnetite (Fe 3 O 4 ) or greigite (Fe 3 S 4 ). Typically, there are a few tens of particles in a magnetosome, and with diameters in the range 30-120 nm.
In this size range, each particle is a single magnetic domain, which maximises the magnetic moment per particle. When these are aligned in a magnetosome, the net dipole moment is large. From the point of view of a nanotechnologist, magnetosomes possess some attractive characteristics. Firstly, the nanoparticles are practically monodisperse, which is quite difficult to achieve in vitro. Biomineralisation is well known to produce inorganic structures with exquisite control. Secondly, the net magnetic dipole moment on a magnetosome is large. This is difficult to achieve with single particles, because although the magnetic moment scales with the volume of a particle, there is a tendency for multiple, disoriented magnetic domains to form in large particles. Alternative strategies such as forming multicore magnetic particles have been devised [4], but still, the dispersion and orientation of smaller particles in a nonmagnetic matrix introduce significant complications. Finally, genetic modifications could provide a means of 'tuning' natural organisms, and harnessing the control and efficiency of biomineralisation. The ability to produce magnetic nanoparticles with well-controlled sizes and shapes, and in particular as aggregates with large net magnetic moments, could lead to optimum materials for certain technological applications.
In medical applications, magnetic particles can be used for the detection and treatment of diseased tissue [5,6]. The accumulation of injected, functionalised magnetic particles in fibrous tissue, typical of tumours, can be detected with pickup coils, meaning that the subsequent surgery can be targeted, and with minimal loss of healthy tissue. In hyperthermia, the application of an AC magnetic field of the right frequency causes dissipative heating of the particles [7], which triggers cell death in the region where the particles have been injected. In this example, the key material property governing the heating rate is the frequency-dependent DMS χ(ω). In response to a weak magnetic field H(t) = H 0 e −iωt with angular frequency ω, the magnetisation is where χ(ω) = χ ′ (ω) + iχ ′ ′ (ω) has real (in-phase) and imaginary (out-of-phase) components. A 'weak' magnetic field here means that the magnetisation is within the linear-response regime. In the classic Debye theory [8,9], in which interparticle interactions are ignored, the DMS is given by where χ(0) is the static magnetic susceptibility, and τ rot is the Brownian rotation time of the particle. (The Debye theory of dielectric response is completely analogous to that for magnetic response.) The heating rate is proportional to χ ′ ′ (ω) [7], and this function has a single maximum at ω = τ −1 rot . For an isolated spherical particle of diameter σ, in a liquid with viscosity η at temperature T, the Brownian rotation time is where k B is Boltzmann's constant. This particular expression applies only to an isolated spherical particle. Interparticle interactions change the functional form of χ(ω), and increase the effective rotation times due to chain-like correlations between the particles, so that τ rot > τ B [10][11][12].
The dynamic magnetic response of suspensions of magnetosomes has been measured in experiments, and the application to hyperthermia treatments has already been considered at length [13][14][15][16][17][18][19][20][21][22][23][24][25][26]. Importantly, the heating efficiency of magnetosomes can be significantly greater than that of normal colloidal magnetic nanoparticles [27,28]. One of the reasons for the enhanced heating effect is that the dipole moments on the constituent magnetic nanoparticles are aligned, giving a large net dipole moment, and hence a large magnetisation response to an external magnetic field.
From a modelling perspective, the structural and static magnetic properties of chains of magnetic nanoparticles have been studied extensively using computer simulations, notably by Kantorovich et al [29][30][31][32][33][34][35][36][37]. In these investigations, not only linear chains of magnetic particles have been considered, but also more exotic architectures such as Y-shaped, X-shaped, and ring SMPs. These studies also address the clustering of SMPs, and how that is affected by architecture, external magnetic field, etc.
Turning to the dynamic magnetic properties of magnetosomes and similar structures, Belovs and Cēbers have predicted some interesting features in χ(ω) [38]. Using a model of a semiflexible dipolar filament, it was found that the bending motions of the filament (and hence the relaxation of the magnetisation) driven by thermal fluctuations at equilibrium lead to a high-frequency decay χ(ω) ∼ ω −3/4 . The theory has to be supplemented with a description of the low-frequency rotation of the filament as a whole, to give a complete prediction of χ(ω). Note that the motion is Brownian, and that a full treatment of the hydrodynamic interactions mediated by the suspending medium would require a more detailed theoretical treatment, numerical schemes, etc. The predicted high-frequency behaviour was tested in experiments on magnetosomes, and the results were consistent with the theory, at least over two decades in frequency [39]. The nonlinear dynamics of magnetic filaments in AC magnetic fields have been shown to lead to a wide range of unusual phenomena [40] including alignment perpendicular to the field [41], oscillating S-, U-, and Z-shaped configurations [41,42], fourarched configurations [42], buckling [43], and microswimming behaviour [44].
There has not yet been a focused study on the DMS of magnetosomes with computer simulations. The aim of the current study is to fill in this gap. The DMS of noninteracting magnetosomes is computed using BD simulations. The specific focus is on how the DMS varies with the number of particles N within the magnetosome. As N is increased, it is anticipated that the DMS changes from the well-known Debye-like function towards something resembling that predicted by Belovs and Cēbers [38]. In addition, the characteristic rotation time should increase with increasing N, according to a range of theoretical predictions [45][46][47][48].
The rest of the article is arranged as follows. The model and simulation methods are described in section 2. The results are presented in section 3, and are organised according to the radius of gyration (section 3.1), the static magnetic susceptibility (section 3.2), and the DMS (section 3.3). Section 4 concludes the article.

Model and simulation methods
A schematic diagram of the magnetosome model is shown in figure 1. It consists of a chain of dipolar soft spheres, each with diameter σ, linked by FENE bonds, and with a bond-bending potential. Each sphere carries a point dipole µ i at its centre.
A soft-sphere interaction operates between all particles, and is given by the purely repulsive WCA potential [49] where r is the centre-centre distance, and ϵ sets the energy scale. The WCA potential is equivalent to the LJ potential cut and shifted at the minimum at r = 2 1/6 σ. The dipole-dipole interaction between two particles i and j is given by [50] where µ 0 is the vacuum permeability, and r = r j − r i is the centre-centre separation vector between the particles. This is the effective interaction between two homogeneously magnetised spheres [50]. Neighbouring spheres in the chain are bonded by FENE bonds, with potentials [51] where K FENE is the spring constant, and the maximum separation is R 0 . When r ≪ R 0 , U FENE ≈ 1 2 K FENE r 2 . Finally, the bond-bending potential is given by the cosine function [52] U bend (θ) = K bend (1 + cos θ) , where K bend is the spring constant, and θ is the angle between neighbouring bonds in the chain. At angles close to the equilib- The bending potential mimics the effects of the membrane encapsulating the magnetic nanoparticles.
In the current work, the temperature is set such that k B T = ϵ, the dipole moment such that µ 0 µ 2 /4πσ 3 ϵ = 4, K FENE σ 2 /ϵ = 30, R 0 /σ = 1.5, and K bend /ϵ = 50. The temperature scale is not too important by itself, as it mainly controls the extent to which the soft spheres can overlap. The FENE and bendingpotential parameters are conventional for coarse-grained simulations of semiflexible polymers [51]. An important physical parameter is the dipolar coupling constant between the particles, which is set here to λ = µ 0 µ 2 /4πσ 3 k B T = 4. This is around the point where unbonded spheres at low concentration spontaneously form chains [53]. This means that the dipoles on bonded spheres will remain aligned, and will not undergo reversal within the frame of the chain. Even a reversal of two Schematic diagram of the magnetosome model. The particle diameter is σ, the dipole moment on particle i is µ i , FENE bonds link neighbouring particles in the chain, and θ is the angle between neighbouring bonds in the chain. spins would require a prohibitive amount of energy: if the energy of the configuration →→ is approximately −2λk B T, and the energy of the 'transition-state' configuration ↑↑ is approximately +λk B T, then the flip rate is proportional to e −3λ . With λ = 4, the activation energy is around 12k B T, and hence the reorientation mechanism for the net dipole moment on the chain is by Brownian rotation of the chain as a whole. A few tests were carried out with λ = 9 and 16, and these will be mentioned in section 3.3. Even so, these values of λ are small compared to the possible values in real magnetosomes with large ferromagnetic or ferrimagnetic nanoparticles. From a computational perspective, smaller values of λ give smaller magnetic forces, and hence larger time steps can be used in integrating the equations of motion.
BD simulations were carried out with LAMMPS [54,55]. As described in detail elsewhere [11,12], BD is actually approximated by Langevin dynamics with a high friction coefficient. The Langevin equations of motion for particle i are [56,57] where superscripts T and R denote translation and rotation, respectively. For the translation of the particle, m is the particle mass, v i is the particle velocity, F i is the net interaction force acting on the particle, ζ T is the translational drag coefficient, and ξ T i is a random force mimicking the Brownian forces arising from the solvent. For the rotation of the particle and its dipole moment, I is the inertia tensor, ω i is the angular velocity, T i is the net interaction torque on the dipole, ζ R is the rotational drag coefficient, and ξ R i is a random torque. Concerning the white-noise, random forces, they obey the statistics ⟨ξ i ⟩ = 0 and ⟨ξ i (t) · ξ j (t ′ )⟩ = 6k B Tζδ ij δ(t − t ′ ), where superscripts T and R are implied. In LAMMPS, each particle is treated as a homogeneous solid sphere, so that I = mσ 2 /10, and the random forces and torques are treated consistently to give the required relationship between the translational and rotational diffusion coefficients of an isolated spherical particle, as expected from the Stokes-Einstein(-Debye) relations: All simulations were carried out using LJ units, denoted with an asterisk ( * ). In LAMMPS, the drag is controlled by a damping time, which is set to damp = 1/20. As shown in [11,12], this gives overdamped, effectively BD with a rotation time τ * B = 1/[6T * × damp] = 10/3 [11,12]. The Langevin equations of motion were integrated with a time step δt * = 0.005. Each simulation was carried out in a cubic box with sides at least twice as large as a perfectly aligned chain, and of course, no long-range interactions were required since the study is restricted to noninteracting magnetosomes. The length of the simulations ranged from 1.67 × 10 6 time steps (or 2.5 × 10 3 τ B ) for N = 1, up to 1.00 × 10 9 (1.50 × 10 6 τ B ) for 16 ⩽ N ⩽ 64, after full equilibration. The effective Brownian relaxation time of the chain, τ rot , increases with increasing N, and so longer runs were required for longer chains. As per theoretical predictions [38], hydrodynamic interactions which would be mediated by the suspending liquid are ignored.
The DMS was obtained through the linear-response route by first calculating the magnetisation autocorrelation function where is the instantaneous magnetisation of the chain, and the angled brackets ⟨. . .⟩ denote an ensemble average. The dynamic susceptibility is given by where χ(0) is the static susceptibility. In the case of noninteracting magnetosomes, χ(0) would be given by where n is the number of magnetosomes per unit volume, but this is arbitrary because a noninteracting, or isolated, magnetosome is being considered. What can be calculated, and will be presented below, is the ratio of the static susceptibility to what it would be if the N particles were not bonded to one another, not interacting, and contained in the same volume. This latter quantity can be called the 'ideal' susceptibility, and the ratio of χ(0) to χ id (0), is given by If the orientations of the N dipoles in a chain were completely uncorrelated from one another, then ⟨M 2 ⟩ = Nµ 2 , and χ(0)/χ id (0) = 1.  The points are from BD simulations, the black line is a fit using (17), and the red dotted line is at Rg = 0. The statistical uncertainties in the BD data are smaller than the symbol size.

Results
Snapshots from simulations of magnetosomes with N = 8, 16, 32, and 64 are shown in figure 2. These images illustrate that the combination of the bending potential and the dipole-dipole interactions renders the chains semiflexible.

Radius of gyration
The radius of gyration of a magnetosome is shown as a function of N in figure 3. Instantaneous values were calculated using the formula and then averaged over the simulation. A suitable fitting function was found to be where M = 3.017 (41), and ν = 0.9026 (45). Hence, R g scales not-quite-linearly with N, reflecting the semiflexibility of the chain.  (18), the dashed black line is a fit using (19), and the red dotted line is at χ(0)/χ id (0) = 1. The statistical uncertainties in the BD data are smaller than the symbol size.

Static magnetic susceptibility
The static magnetic susceptibility compared to the ideal value (15) is shown as a function of N in figure 4. Note that with N = 1, χ(0)/χ id (0) = 1. The BD simulation data are fitted well with a theoretical expression for magnetic filaments [58,59], where κ is a parameter that characterises the rigidity of the chain. When κ = 0, χ(0)/χ id (0) = 1, corresponding to the N dipole orientations being completely uncorrelated with one another. When κ = 1, χ(0)/χ id (0) = N, which corresponds to perfect alignment of the dipoles at all times (|M| = Nµ). Fitting (18) to the BD simulation data gives κ = 0.9718 (11). This value being close to 1 confirms the strong alignment of the dipoles within the chain, as found in real magnetosomes. Another way of analysing the simulation data is to assume that, since the dipoles are aligned in a magnetosome, ⟨M 2 ⟩ is proportional to the square of the radius of gyration. Hence, χ(0)/χ id (0) was also fitted with the equation where P = 0.771 (36), and ν = 0.886 (10). The fit is shown in figure 4. The value of ν is consistent with the value obtained from R g , which confirms that the basic reasoning is correct, although (18)

Figures 5(a) and (b)
show, respectively, the real and imaginary parts of the DMS. All of the data are shown scaled by the static susceptibility, χ(0), as per (13). Each curve is labelled with the number of particles N in the magnetosome. Note that the data do show some statistical scatter, particularly with large values of N. With large values of N (N ⩾ 32), C(t) decays extremely slowly, and computing the long-time decay accurately is challenging. In carrying out the Fourier transform of C(t) (13), a Blackman windowing function was used in order to minimise low-frequency truncation errors; this should not affect the high-frequency portion very much. The variations in peak heights in χ ′ ′ (ω) in figure 5(b) are probably not meaningful. Figure 5(a) shows that the drop in χ ′ (ω) occurs at lower frequency with larger N. With N = 1, τ rot = τ B , and hence χ ′ (ω) = 1 2 at ωτ B = 1; see (2) and (3). As N is increased, the curveshifts to lower frequencies, signalling a growth in the effective Brownian rotation time. Figure 5(b) shows χ ′ ′ (ω) for different values of N. With N = 1, the peak χ ′ ′ (ω) = 1 2 is at ωτ B = 1, and as N is increased, the maximum shifts to lower frequencies. It is not immediately obvious from these plots that  (2), and the Belovs-Cēbers theory (20). The dimensionless frequency is ωτ bend . For the purposes of comparison, the Debye expression is shown with τrot = 13τ bend /55 440, so that the low-frequency dependence is the same as the Belovs-Cēbers prediction. the shape of the DMS changes with increasing N, but the data can be used to test theoretical predictions.
Belovs and Cēbers derived the following expression for the DMS of a semiflexible, dipolar filament [38,39].
Here, τ bend is the time scale for thermal fluctuations in the magnetisation of the filament due to bending motions only, and p is the complex conjugate of the variable used in equation (4) of [39], so that The second term in (20) removes a divergence in χ ′ ′ (ω) at low frequencies, and yields the correct form in that limit [38,39]. If, as p → 0, and Note that p 4 = iωτ bend , and hence with the diverging p −4 term removed, and multiplication by 420, the p 4 term means that χ ′ ′ (ω)/χ(0) ≈ 13ωτ bend /55 440 at low frequency.
It is useful to plot (2) and (20) together, to highlight the differences. Figure 6 shows this comparison, with the curves plotted as functions of ωτ bend . The peak in the Belovs-Cēbers function is at ωτ bend ≃ 3877. If one identifies the peak position with ωτ rot = 1, then τ bend /τ rot = 3877. For the purposes of comparison, the Debye curve is shown with τ rot = 13τ bend /55 440, so that the low-frequency dependence is the same. These plots confirm that, at low frequency, the results are the same by construction, while at high frequency, the Belovs-Cēbers expression peaks a little earlier, and has a broader tail extending to higher frequencies. One of the key predictions from the theory is that χ(ω) ∼ p −3 ∼ ω −3/4 at high frequencies.
To summarise the three relevant time scales, τ B is the Brownian rotation time for an isolated, spherical particle (3). τ rot is the effective rotation time of a magnetosome, which can be identified from the peak position in χ ′ ′ (ω), and it is larger than τ B because a magnetosome is larger than a sphere. τ bend is the typical time scale for fluctuations in the magnetisation of a magnetic filament due the bending, and it is much larger than τ rot because net rotation of the magnetosome leads to faster reorientation of the magnetisation than bending motions alone. Figure 7(a) shows log-log plots of χ ′ ′ (ω) from BD simulations with all values of N, along with the power law χ ′ ′ (ω) ∼ ω −3/4 . As N is increased, the high-frequency portions begin to show the predicted power-law behaviour, and with N = 32, 48, and 64, this occurs over several decades in frequency. With these large values of N, the simulation results at very high frequency (ωτ B ⩾ 10 −1 ) are noisy because of the sampling frequency when calculating the magnetisation autocorrelation function C(t). As noted above, the windowing function does not strongly affect the high frequency portion (10 −4 ⩽ ωτ B ⩽ 10 −1 ). The Belovs-Cēbers prediction of the behaviour at high frequencies also applies to the real part of χ(ω), but that is practically impossible to study with the BD simulation data,  Fitting parameters for magnetosomes with N = 64, and λ = 4, 9, and 16. τ B is the Brownian rotation time of an isolated spherical particle, τrot is the Brownian rotation time of a magnetosome, τ bend is the bending time for a magnetic filament, and ω 0 and a are parameters in a power law (see text). The functions used in fitting BD simulation data are given under each parameter. 'Average' indicates the mean value of the parameters, and 'all' means a fit to all of the data simultaneously.  (12) because χ ′ (ω) ≪ χ ′ ′ (ω) (see figure 6), and the signal-tonoise ratio is too poor. Figure 7(b) shows log-log plots of χ ′ ′ (ω) from BD simulations with N = 64, and λ = 4, 9, and 16. The different sets of results are very similar, and from the form of the plots, it is clear that the increase in χ ′ ′ (ω) at low frequency is more rapid than the decrease at high frequency. The high-frequency data follow power-law behaviour over almost three decades in frequency.
The BD simulation results were fitted with the Debye (2) and  expressions in the range ωτ B ⩽ 1 × 10 −1 . Parameters for each value of λ, averages of those parameters, and parameters from fitting all of the data simultaneously are given in table 1. The typical Debye rotation time is τ rot ≃ 2.8 × 10 4 τ B , and the typical bending time is τ bend ≃ 1.1 × 10 8 τ B . Plots of the functions with these typical parameters are shown in figure 7(b).
The Debye theory (2) predicts that at low frequency, χ ′ ′ (ω) ∼ ω, and at high frequency, χ ′ ′ (ω) ∼ ω −1 . It is clear that the Debye theory is only correct at low frequency. While the Belovs and Cēbers theory is by no means accurate over the whole frequency range, it does describe the non-trivial highfrequency scaling correctly. Fitting a function (ω 0 /ω) a to the high-frequency BD simulation data over the range 5 × 10 −4 ⩽ ωτ B ⩽ 1 × 10 −1 gives the parameters listed in table 1. From the fits for each value of λ, the exponent a is in the range 0.7-0.8. Averaging the individual exponents gives a = 0.758 (12), and fitting all of the data simultaneously gives a = 0.755 (12). These values of a are consistent with the Belovs and Cēbers prediction of a = 3 4 , and so the Debye theory in this regime can be ruled out.
A characteristic rotation time, τ rot , of a magnetosome can be extracted from the position of the maximum in χ ′ ′ (ω), at which ωτ rot = 1. (Another option is the time it takes for C(t)/C(0) to reach 1/e, but purely exponential decay cannot be assumed at all times, at least with long magnetosomes.) The maximum was found by fitting a third-order polynomial in the region of the peak. The results with λ = 4 are shown as a function of N in figure 8. There have been many theoretical predictions for how the rotation time of a cylindrical particle should increase with increasing aspect ratio [45][46][47][48]. The aspect ratio is L/D, where L is the length and D is the diameter of the particle. Roughly speaking, if D ≃ σ, and L ≃ Nσ, then L/D ≃ N. The results in figure 8 are compared to an expression similar to that used by Yan et al [48], and with no fitting parameters. The comparison with the BD simulation data is excellent. Finally, the Belovs-Cēbers theory predicts a maximum in χ ′ ′ (ω) at ωτ bend ≃ 3877, and hence τ bend /τ rot = 3877. This ratio, from the Belovs-Cēbers and polynomial fits to the BD simulation data, is given in table 1. The ratios for different values of λ are about 20% higher than the theoretical prediction. Matching the low-frequency behaviours of the Belovs-Cēbers and Debye expressions gives τ bend /τ rot = 55 440/13 ≃ 4265. The deviation of the simulation results from this value is about +9%.

Conclusions
The properties of noninteracting, model magnetosomes have been calculated using BD simulations. Each magnetosome contains N strongly interacting ferromagnetic particles, connected in a chain with nonlinear elastic bonds, and with a bond-bending potential. The radius of gyration, static magnetic susceptibility, and DMS have been computed for N = 1-64. The radius of gyration follows a scaling law R g ∼ N 0.9 , and the scaling of the static magnetic susceptibility with N corresponds to a stiffness parameter κ ≃ 0.97. With increasing N, the DMS crosses over from the normal Debye form to a more complex function. Focusing on the imaginary (out-ofphase) part of the susceptibility, while the Debye prediction is that χ ′ ′ (ω) ∼ ω −1 at high frequency, the Belovs-Cēbers prediction is that χ ′ ′ (ω) ∼ ω −3/4 . At low frequency, both theories predict that χ ′ ′ (ω) ∼ ω. The high-frequency simulation data are much closer to the Belovs-Cēbers scaling law when N is large (N ⩾ 32). This reflects the fact that a long chain much more closely resembles the homogeneously magnetised, semiflexible dipolar filament considered by Belovs and Cēbers. The peak in χ ′ ′ (ω) shifts to lower frequency as N is increased, reflecting a growing Brownian rotation time, τ rot , for the magnetosome as a whole. It was confirmed that the dependence of τ rot on N follows well-established theoretical results.
The simulations were carried out in the absence of hydrodynamic interactions, i.e. the Brownian forces exerted by the suspending liquid on the constituent particles were uncorrelated with one another. This was appropriate for the purpose of comparing simulations with theory. In reality, the motion of one part of the magnetosome causes a disturbance in the velocity field of the suspending liquid, which affects, and is hence correlated with, the force exerted by the medium on another part of the magnetosome. Techniques such as the lattice-Boltzmann method could be used to assess the effects of such hydrodynamic interactions. Other avenues for future research include magnetosomes with different kinds of engineered structures such as X-or Y-shaped nanoparticle aggregates, the effects of interactions between different magnetosomes, and the effects of crowding and confinement that occur in cells and fibrous tissue.
The practical relevance of such results has already been recognised. In magnetic hyperthermia therapy, an AC magnetic field with angular frequency ω causes heating in a magnetic medium at a rate proportional to χ ′ ′ (ω), and localised heating can be used to treat diseased tissue. The present simulation results confirm how χ(0) and χ ′ ′ (ω) increase with N, and how the optimum frequency for heating decreases with N.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).