Optimized shapes of magnetic arrays for drug targeting applications

Arrays of permanent magnet elements have been utilized as light-weight, inexpensive sources for applying external magnetic fields in magnetic drug targeting applications, but they are extremely limited in the range of depths over which they can apply useful magnetic forces. In this paper, designs for optimized magnet arrays are presented, which were generated using an optimization routine to maximize the magnetic force available from an arbitrary arrangement of magnetized elements, depending on a set of design parameters including the depth of targeting (up to 50 mm from the magnet) and direction of force required. A method for assembling arrays in practice is considered, quantifying the difficulty of assembly and suggesting a means for easing this difficulty without a significant compromise to the applied field or force. Finite element simulations of in vitro magnetic retention experiments were run to demonstrate the capability of a subset of arrays to retain magnetic microparticles against flow. The results suggest that, depending on the choice of array, a useful proportion of particles (more than 10%) could be retained at flow velocities up to 100 mm s−1 or to depths as far as 50 mm from the magnet. Finally, the optimization routine was used to generate a design for a Halbach array optimized to deliver magnetic force to a depth of 50 mm inside the brain.


Introduction
Magnetic drug targeting (MDT) has recently become a topic of interest among researchers due to its potential to localize and retain therapeutic agents efficiently in a target region, which has possible applications for the treatment of a range of diseases including cancer [1][2][3][4][5] and damaged blood vessels [6][7][8]. MDT using permanent magnets is advantageous because static magnetic fields and forces can be applied inside the body without being attenuated by tissue or posing a risk of magnetic hyperthermia [9,10]. There are, however a number of challenges to overcome before the technique can be considered clinically viable [9,11]. A major issue is that magnetic fields and, hence, forces decay rapidly with distance, limiting applications to relatively shallow targets in the human body [12][13][14]. Additionally, the applied magnetic force must overcome the hydrodynamic drag force of blood before a useful quantity of agent can be captured and retained against the flow of the circulatory system [15][16][17][18].
MDT delivery systems usually consist of a therapeutic agent contained within a bio-compatible carrier functionalized or loaded with superparamagnetic iron oxide nanoparticles, and much work has focused on increasing the magnetic moment of these carriers [19][20][21][22]. Mesoscopic magnetic carriers (nanometres to microns) are particularly interesting because of a favourable ratio between magnetic and Stokes' drag forces allowing for improved accumulation [10], and the ability to tailor multi-modal composite carriers that may and are responsive to external stimuli for imaging [26] and controlled release [27,28]. However, while it is well understood that the carrier formulation needs to be optimized for the application [5,11,29,30], there is increasing recognition that the external magnet system also needs to be tailored to the requirements and constraints of a given application,  accounting for the depth and physiological flow regime at the target [14,17,[31][32][33][34][35][36].
In our previous study [37], we developed an optim ization routine to determine the configuration of an assembly of inexpensive, readily available cubic permanent magnet elements offering the maximum field gradient at a given distance from the array. The aim of the present study was to significantly expand this approach to the design of magnetic arrays consisting of elements of arbitrary shape. We show that the resultant arrays are capable of generating almost two to three times as much magnetic force as arrays constructed using cubic elements for the same volume of magnetic material, depending on the optimization distance.
The difficulty of assembling arrays consisting of multiple permanent magnet segments due to the repulsive dipole forces that arise in some configurations is also considered and addressed. Designs of shapes generated using a uniform magnetization are proposed as appropriate for the soft ferromagnetic core of 'open-pot' electromagnets, such as that reported by Alexiou et al [38]. Finite element simulations of a subset of arrays are performed to demonstrate capture of magnetic particles in a range of physiologically relevant flow velocities and at depths up to 50 mm from the magnet surface. Finally, the versatility of the optimization routine is demonstrated in the form of a design of a Halbach array specifically tailored to actuate and retain magnetic particles against flow at different tissue depths inside the brain.

Model of magnetic force
A general expression for the magnetic force, F, on a single domain superparamagnetic particle with a moment of V M B ( ) µ = is given by where M is the magnetization of the particle, which depends on the field, V is the volume of the particle and B H 0 µ = is the magnetic flux density, proportional to the applied field, H. As the particle is superparamagnetic, it is assumed that M and B are parallel. The magnetization of a superparamagnetic particle can be described using a Langevin function, L y y y coth where M s is the saturation magnetization of the particle, H is the applied field inside the particle and k T B is the product of the Boltzmann constant and the temperature [39][40][41].
The field emitted by an array consisting of an arbitrary configuration of magnetic elements was calculated by breaking the magnet into a 3-dimensional arrangement of evenly distributed point moments, following a method described previously [37]. Each moment emits a dipole field described by is the point moment, M is the magnetization of the permanent magnet, V d is the volume occupied by the point and ′ r is the position vector relative to the point moment. In the optimization routine described below, the normalized magnetic force due to the field emitted by an array of magnets on a superparamagnetic particle at a position of interest (POI) was calculated. The normalized magnetic force (or force per moment) is given by and has units of T m −1 [38]. When the particle is saturated (M M s = ), the normalized force is equivalent to the field gradient emitted by the array. The superparamagnetic particle considered here has the same saturation magnetization as Fe 3 O 4 at room temperature (M 4.7 10 s 5 = × A m −1 ) and a diameter of 10 nm.
The model was implemented using console applications written in the C programming language (Microsoft Corporation, Redmond, WA, USA).

Optimization routine
An optimization routine was developed to generate designs of arbitrarily-shaped magnet arrays to deliver the maximal normalized force on a particle at the POI (r POI ) given a series of design parameters, including the volume to be optimized, the nominal direction of normalized force (F nom ), the volume of the magnet (V mag ), and the list of allowable magnetization directions contained within the array ( figure 1). An initial array is constructed to occupy the volume to be optimized consisting of both magnetized and non-magnetized elements, with magnetized elements occupying the positions closest to the POI. The total volume of the magnetized elements is limited to V mag at each step using a subroutine described below.
The main routine then starts at the element closest to the POI and tests each allowable magnetization orientation, retaining the one that results in the best value of the optimized parameter, generated by the whole array at the POI. The process is then repeated for the next closest element until all elements in the array have been treated. At this point, convergence is tested by comparing the attained array to the configuration of the starting array. If the routine has changed the array and resulted in an improvement in the optimized parameter, the process is rerun using the attained array as the new starting array and again starting from the element closest to the POI until all elements have been treated. If the routine does not change the array after treating all elements and the optimized parameter cannot be improved, the array is considered optimized.  Whenever the combined volume of all elements with a non-zero magnetization exceeds the V mag parameter, a subroutine is performed in order to find and demagnetize the element that makes the least contribution to the normalized force. As the force depends on the gradient of the total field generated by the array at the POI, it cannot be assumed that this element is the element furthest from the POI. To find the element to demagnetize, each magnetized element is temporarily replaced by a non-magnetized element of the same volume and for the remaining array is recorded. The element that makes the least difference to the optimized parameter when replaced by a non-magnetic element is demagnetized.
A collection of optimized arrays were generated by varying different parameters in the optimization routine, including the volume of the magnet, the distance between the magnet and the POI, the direction of force and the set of allowable magnetizations. Unless otherwise specified, the initial magnet array was constrained to a 0.2 0.2 0.1 × × m 3 optimization volume positioned directly below the x-y plane, with the POI set along the z-axis above the x-y plane. The default element density was 4 cc −1 . Three sets of allowable magnetization directions were investigated, within which all possible magnetization vectors had the same magnitude as that of an N52 grade NdFeB permanent magnet (1.14 10 6 × A m −1 ). The first set (uniform magnetizations) contained a magnetization vector aligned with the z-axis and a zero vector (totaling two possible configurations). The second set (orthogonal magnetizations) contained six vectors pointing in the positive and negative of each orthogonal direction, along with a zero vector (seven possible configurations). The third set (diagonal magnetizations) contained all vectors in the orthogonal set, along with all possible corner and edge diagonal directions and a zero vector, totaling 27 possible magnetization configurations.
The routine returns the position and magnetization of all magnetized elements at optimization. Figure 2 shows how a resultant arrangement of magnetization vectors can be interpreted and converted into a design of constructable shapes and dimensions, particularly when the output is approximately cylindrically symmetrical (which is often the case). This is done by merging regions with the same magnetization into individual segments. The difficulty of assembling large permanent magnet segments that, in some configurations, can be strongly repulsive, must be considered [37] and is quantified here in terms of the internal magnetic potential energy, U int , which is calculated by summing B i int µ − ⋅ for each element,  where B int is the field generated at the dipole position by the array without the merged magnet segment containing the element. Larger values of U int are interpreted as arrays that are more difficult to assemble.

Finite element simulations
Finite element modeling was performed using COMSOL Multiphysics 5.0 (COMSOL, Inc, Burlington, MA, USA) to assess the suitability of a subset of optimized arrays for magnetic drug targeting applications via particle tracing simulations. Particles with the same properties as magnetically loaded polystyrene microbeads were simulated in a laminar flow within a straight channel primed with water flowing above the magnet in the x-direction. Capture efficiency was determined by quantifying the proportion of particles that accumulated in the region above the magnet after a simulation time of 200 s. The field emitted by a given array was calculated using the 'magnetostatics, no currents' interface of the 'AC/DC module' after constructing the geometry of the array in a 3D model following the method described for figure 2. The capture efficiency and microbead accumulation was then approximated to first-order [42] in a simplified 2D geometry by first modeling the flow profile in a 3 mm wide channel using the 'laminar flow' interface of the 'CFD module' (setting a no-slip boundary condition at the wall and zero outlet pressure), and then using the 'particle tracing module' to solve the trajectories of particles in flow under the influence of a drag force, gravity and a magnetophoretic force described by the following three equations respectively: where the various parameters are given in table 1. The field in (5c) was taken as an interpolated function of the solution to the 3D array model in the x-z plane (figure 3). The magnetic permeability of microbeads, r,p µ was set to 1 χ is the magnetic susceptibility, and M H ( ) is described using (2) and assuming an effective superparamagnetic cluster diameter of 10 nm to account for the fact that particles approach magnetic saturation when exposed to the particularly intense fields emitted at the face of an optimized array. The saturation magnetization of microbeads was set to , where α is the volumetric ratio of superparamagnetic Fe 3 O 4 in microbeads. The density of particles was given by 1

Volume-dependent optimizations
Optimizations were performed to generate magnet designs of different volumes between 10 and 1000 cc, using one of the three magnetization vector sets described in section 2.2, with the POI set at 20 mm along the z-axis and F nom directed towards the magnet (in the negative z direction). Figure 4 shows 2D element maps of a subset of optimized designs with V mag constrained to 50 or 500 cc. Using a larger magnetization vector set results in arrays that are more tightly confined closer to the POI. Typically, uniform arrangements have tapered tips on the 'front' face of the magnet (closer to the POI), causing the field to decay more rapidly in the region close to the POI (increasing the B ( ) ∇ component of (4)). This isn't the case for the overall Halbach arrangements (consisting of orthogonal and diagonal magnetization sets), but the central magnet segment (   on position along the z-axis for a set of arrays with magnet volumes of 10 and 1000 cc. In each case, past a certain distance (usually about 15 mm) the field and force profiles decay approximately exponentially with distance (straight on a lin-log graph), highlighting the difficulty of applying useful magnetic forces over a long spatial range. The diagonal arrangements are able to exert significant normalized forces very close to the face of the magnet, exceeding 100 T m −1 , even for the 10 cc magnet, but the uniform shapes tend to perform better at long range, with the 10 and 1000 cc uniform magnet applying greater fields at z = 50 mm than Halbach arrangements of the same volume and not decaying as quickly in normalized force. Over the entire displayed range, the diagonal Halbach arrangements are superior to orthogonal designs with the same volume.
The effect of changing magnet volume on the field and force generated at the POI is displayed in figure 6. This indicates that for each magnetization set, B z POI ( ) and F z MV POI s ( )/ increase approximately logarithmically with V mag and, for designs with orthogonal magnetizations, a factor of ∼5 increase in volume is required to raise the normalized force by 10 T m −1 . Notably, the 100 cc orthogonal array produces a normalized force of 20.3 T m −1 at the POI, almost twice the force of a double layer cubic element array at the same distance (11.9 T m −1 ) reported in reference [37], which was optimized using the same set of parameters.
The internal magnetic potential energy, shown in figure 7(a) was calculated following the method described in section 2.2 as a metric to gauge the difficulty of assembling an optimized array from individual segments of uniform magnetization. Uniformly magnetized shapes have a U int of 0 J assuming they are machined from a single piece and not assembled from smaller elements that are repulsive in certain configurations. As the number of possible magnetization vectors increases and, thus, the number of segments required to assemble an array, the internal magnetic potential energy also rises. This demonstrates the main disadvantage of designs using diagonal magnetizations; while arrays utilizing diagonal configurations tend to result in the most intense field and force values, their assembly is complicated by the fact that many of the configurations between neighbouring segments are repulsive. The potential energy density in the 1000 cc array with diagonal magnetization vectors is 3.3 10 5 × J m −3 . The required computation time to execute the optimization routine for each design on a computer with an Intel(R) Core(TM) i7 processor and 8 GB RAM is shown in figure 7(b). Uniformly and orthogonally magnetized arrays with V 100 mag ⩽ cc and an element density of 4 cc −1 can be optimized in under an hour. Diagonally magnetized arrays require significantly more time to optimize.

Position of interest dependent optimizations
A set of optimizations was performed with the magnet volume constrained to 100 cc and the nominal direction of force fixed towards the magnet (pull force). The position of interest was varied along the z-axis and designs were generated for each of three possible magnetization sets described previously. Figure 8 shows the resultant designs with the POI set at distances of 5 and 50 mm away from the magnet. For each magnetization set, increasing z POI yields designs that are more tightly confined to the x-y plane, resulting in flatter, disc shaped volumes. When the POI is very close to the face of the The relatively high proportion of elements off-axis that are magnetized away from the POI results in extremely high field and force values very close to the face of the arrays, as shown in the profiles of B and F M V s / in figure 9. Both types of Halbach arrays optimized for z 5 opt = mm are able to obtain a field of 1.6 T and a field gradient in excess of 300 T m −1 at short range, and are even capable of applying normalized forces greater than 100 T m −1 as far as 7 mm away. These field gradients are remarkable, and are more than twice as forceful as unresolved (i.e. difficult to assemble) cubic, pull arrays of the same volume reported previously (maximum B ( ) ∇ of 139 T m −1 ) [37], and almost three times as forceful as the corresp onding resolved arrays which were determined to be easier to construct (capable of applying 124 T m −1 ). To our knowledge, the only other magnetic systems capable of applying field gradients of several hundred T/m over millimetre length scales that have been considered for MDT are based on superconducting magnets [43,44]. However, the field profiles in figure 9(a) are not persistent and drop off fairly rapidly, with the magnitude of B exhibiting minima at about 25 mm, coinciding with a change in direction of field from positively aligned with the z-axis to negatively aligned. This results in very small push forces from the magnets in the region between ∼25-40 mm, between the two minima in F M V s / ( figure 9(b)). The fields and forces emitted by uniformly magnetized arrangements decay less rapidly as a function of distance than those emitted by orthogonal and diagonal Halbach arrays optimized for the same POI. Similarly, arrays optimized for a further POI also perform better over a longer range of distances than arrays optimized for z 5 opt = mm. The diagonal array optimized for 50 mm is capable of delivering a normalized force greater than 10 T m −1 at 28.5 mm, while the field remains above 0.1 T up to 43 mm away. Figure 10 shows how the field and normalized force vary for different arrays at the POI, as a function of POI. The behaviour of B z POI ( ) highlights how the advantages of using Halbach arrays over uniform magnets diminish at long ranges, particularly for orthogonal arrangements at z 50 POI = mm, where the improvement over the uniform design is indiscernible. Figure 10(b) shows that, while diagonal arrangements apply superior forces at all POIs, the improvement over orthogonal arrays is relatively small, particularly when  the POI is set close to the magnet. When z 5 POI = mm, the performance of the diagonal arrangement is 10% better than the orthogonal array, while, for z 50 POI = mm arrays, using a diagonal magnetization set results in a 20% greater force than the orthogonal design.
The behaviour of the field and force in a x-y and x-z plane from an orthogonal magnet design with z 20 POI = mm is displayed in figure 11. It is noted that, while the force in the x-y plane 5 mm above the surface of the magnet is directed towards its centre axis, the strongest forces coincide with the regions where the magnetization changes on the upper surface of the design.

Direction of force dependent optimizations
The optimization routine was used to investigate how optimized designs vary with different nominal directions of force. V mag was fixed to 100 cc and the POI was set at 20 mm along the z-axis, while the angle between the nominal direction of force and the negative z-axis, labeled θ, was rotated through the x-z plane (this convention was chosen so that 0 θ = results in the nominal direction of force pointing toward the magnet, optimizing for a pull force, and 180 θ = coincides with a nominal direction of force away from the magnet, to maximize the push force at the POI).
A subset of resultant designs consisting of uniform and orthogonal arrays is displayed in figure 12. For the uniform array with 90 θ = , in order to obtain a component of force in the −x-direction along the z-axis, the array splits into two parts, with most of the volume of the magnet occupying the −x, −z quadrant. When 180 θ = , a toroid shape results, centred around the z-axis. With this geometry, a push force results at the POI because the field close to the face of the toroid is negative along the z-axis (the cyan line in figure 13(a)). This field changes direction at a distance along the z-axis dictated by the geometry of the toroid; the transition from negative to positive field results in a push force away from the local minimum in field ( figure 13(b)).
Setting 180 θ = with orthogonal magnetization vectors essentially 'inverses' the pattern in the centre of the array, with a segment magnetized anti-parallel to z occupying the central axis, and surrounding elements magnetized to redirect flux away from this segment and into a toroid with the same magnetization and comparable dimensions to the shape generated for the uniform case.
Close to the face of the array, the field for arrays optimized with 90 θ = is almost parallel to the x-axis, while the field for arrays with θ set to 180 points in the negative zdirection ( figure 13(a)). However, as the distance from the array increases, the field vector rotates to be more closely aligned in the z-direction, either gradually in the case of 90 θ = arrays, or as a sudden transition in the case of push force arrays. Figure 13(b) shows that, even arrays optimized for 180 θ = are only capable of applying a push force over a short distance of ∼15 mm along the main axis, as seen by the range for which the the force points at 180 to the −z-direction (between the two minima in the magnitude of the force profiles). The range of push force is slightly greater for the uniform arrangement, but the orthogonal array gives a slightly greater force. Push force arrays can be useful for noninvasive magnetic injection when physiological flows in the region of interest are low [45,46], but for applications where carriers need to be separated from high blood flow velocities (e.g. partially-occluded and/or injured arteries, or around the leaky vasculature of tumours), high field gradients (more than 10 T m −1 ) that persist over a range of several centimetres are more useful [38,47].
The capability of arrays optimized with different θ values to deliver field or force to the POI 20 mm away is displayed in figure 14. Of interest is the fact that the angle of the relevant vectors at the POI is largely independent of the magnetization set used to generate the design, but diagonal magnetized arrays are consistently able to deliver about twice as much normalized force than the analogous uniform arrangements. However, while the diagonal magnetization set performs the best of the tested sets with all θ, the difference in performance between diagonal and orthogonal arrangements is diminished when a push force is the objective.

Particle tracing simulations
Particle tracing simulations were performed using COMSOL software following the method described in section 2.3, to calculate the trajectories of magnetic microbeads in fluid flow past a subset of optimized magnet arrays. The field profiles for three orthogonal magnet arrays (z 5 POI = , 20 and 50 mm as reported in section 3.2) were calculated by assembling 3D models following the method described in figure 2. Fluid flow velocities and particle trajectories were calculated inside a 2D channel that was 3 mm wide, with the channel centre-line (corresp onding to max fluid velocity) positioned at various distances, z d between 5 and 50 mm above the upper surface of the array. As the laminar flow profiles were calculated in 2D, the velocity along the centre-line of the channel was 1.5 times the nominal, mean flow velocity, which was varied between 2.5 and 250 mm s −1 (corresponding to Reynolds numbers between 8.4 and 840). The lower end of this range is of the same order as blood flow velocities in the cerebral cortex [48,49], while the upper end of this range is comparable with those observed in tumours [50,51].
Mesh independence was determined for the 3D model of the z 50 POI = mm orthogonal array by using the available physics-controlled meshes and increasing the mesh density until no variation was observed in the field profile along the z-axis. Figure 15(a) shows that mesh independence is obtained once the mesh density exceeds ∼2.14 10 6 × elements m −3 and that the finite element calculations agreed well with the dipole model for all meshes. As calculation of magnetic field using COMSOL's interfaces was not computationally intensive on the PC described in section 3.1, the finest available physicscontrolled mesh density was used for all subsequent calculations of field profiles for the other arrays. A similar procedure was followed for the 2D model by calculating the fluid velocity magnitude along the width of the channel ( figure 15(b)), and mesh independence was attained for meshes with a density greater than 9.33 10 7 × elements m −2 . A mesh density of 2.68 10 8 × elements m −2 was used for particle tracing simulations (containing a minimum element size of 2.08 10 6 × − m). Figure 16 shows the capture efficiency of microparticles as a function of flow velocity using different arrays and with the channel set at different distances. In each case, the capture efficiency decays approximately as a power law with flow velocity, typically with an index of about −0.6 (closer to −0.7 when z d = 50 mm). At a channel position of 5 mm, there is very little difference in the total capture efficiency of the z 5 POI = mm and z 20 POI = mm arrays, although there is  a difference in the distribution of captured particles along the length of the channel ( figure 17(b)). At further channel distances, the z 5 POI = mm array is vastly inferior to the other two arrays due to the fact that its applied force decays most rapidly with space, while the array that performs the best is the one that was optimized for that range. Our previous work on magnetic carriers [37,52] has indicated that a capture efficiency of 10% is sufficient to significantly increase the acoustic response detected from retained magnetic microbubbles under ultrasound exposure and this would in turn be expected to generate a therapeutic effect [53]. On this basis an optimized array (i.e. z z d POI = ) would be able to retain a diagnostically and/or therapeutically relevant number of particles in flow velocities of 100 mm s −1 at 5 mm, 25 mm s −1 at 20 mm and 2.5 mm s −1 at 50 mm.
The accumulation distribution was quantified by counting the relative proportion of captured microparticles that were distributed along the length of the channel above the magnet in 5 mm incremental sections. The accumulation distribution inside the channel set 5 mm above the z 20 POI = mm array is     figure 17(a). At all velocities, two peaks emerge, one coinciding with the z-axis, above the centre of the array, and one coinciding with the leading edge of the array, where the field profile has a local shoulder. Figure 17(b) shows how accumulation varies using the different arrays and setting the inlet velocity to 10 mm s −1 . At this range, the z 5 POI = mm array is better able to localize more particles in the region near to the centre of the array, owing to a strong attraction to the particularly intense and narrow peak in field profile in this region along the channel (displayed in the inset). The z 50 POI = mm array exhibits a local minimum in the field around this region, resulting in very few particles being directed to the target.

Consideration of assembly forces
In order to assemble Halbach arrays, the sometimes repulsive dipole forces that arise between neighbouring permanent magnet elements must be overcome. This challenge can be somewhat mitigated by separating neighbouring segments that are in repulsive configurations, but this may result in a compromise in the field and force generated by the array along the z-axis. The designs of the 100 cc orthogonally and diagonally magnetized Halbach arrays optimized for z 50 POI = mm (displayed in figures 8(e) and (f )) were resimulated after introducing a minimum separation distance, d between each segment. Figure 18 shows how introducing a gap up to 10 mm between segments can reduce the internal magnetic potential energy associated with each array, U int which is used to gauge the relative difficulty of assembling a given array. Introducing a gap of 10 mm into the orthogonal array lowers U int by more than 60%, compared with the same design with no gap. Our model allows us to compare this quantity to that of other arrays reported in the literature, such as a Halbach cylinder (108 o.d. 54 i.d. 115 × × mm 3 ) assembled by Cugat et al [54]. Our analysis suggests this design has a U int value of 10.6 J, which is favourably comparable to the values we report for optimized orthogonal arrays with separations greater than ∼6 mm, implying that the challenges associated with assembling these designs can be overcome. The potential energy density of the orthogonal design with d = 10 mm is 6.56 10 4 × J m −3 . The field, B and force, F generated by the orthogonal and diagonal arrays with different separation distances were calculated along the z-axis at positions 5 and 50 mm away from the magnets. Figure 19 shows how B/B 0 and F/F 0 vary for these two arrays at these two positions for different gap sizes, d, where B 0 and F 0 are the field and force respectively when d = 0 mm. Varying d has a much greater effect on the diagonal array than the orthogonal array; a separation of 10 mm only diminishes the field generated by the orthogonal array at 5 mm by ∼11%, and the compromise at 50 mm is even less ( figure 19(a)). Interestingly, the force at 5 mm, displayed in figure 19(b), increases for both arrays for larger d, a consequence of the fact that the field gradient in this region is changing more rapidly due to a faster decay in field. At z = 50 mm, F generated by the diagonal array decays so rapidly with d that, for separations greater than 6 mm, the orthogonal array produces a greater total force than the  diagonal array. This is notable because no other parameter set investigated in the present study yields an optimized arrangement from the orthogonal magnetization set that is capable of applying more force than the analogous design using the diagonal magnetization set.

Example application
The optimization routine can readily be adapted for specific magnetic drug targeting applications by defining parameters related to the position of interest, nominal direction of force and volume of magnetic material, and also the dimensions of the constrained shape of the optimization volume, taking into account the anatomy and physiology (in particular, fluid flow and vessel diameter in the vessel network near the target) of the targeted region. A design volume was set up to optim ize a Halbach array of orthogonally magnetized segments to target and retain magnetic microparticles in the vessel network around the brain at a depth of 50 mm, which has a range of possible applications (for example, magnetic microbubbles can be used for localized opening of the blood-brain barrier to deliver drugs to brain tumours [26]). The helmet-shaped design volume consisted of a hemisphere with an internal radius of 100 mm and a thickness of 15 mm above the x-y plane, and the magnet volume was constrained to 200 cc (equivalent to a magnet weight of 1.5 kg). The direction of force was set at 45 to the z-axis, to provide a component of force that acts against the general direction of blood flow in the region of interest, along with the component of force that pulls towards the magnet.
The resultant design is shown in figure 20. The top view is similar to the designs of orthogonally magnetized arrays discussed in previous sections, skewed to be more weighted towards the −x quadrants to accommodate the diagonal direction of force. The performance of the array along the z-axis is exhibited in figure 21. Past about 25 mm from the inside surface of the magnet, the direction of field and force varies little while the magnitude of these vectors decays approximately exponentially. However, notably the field decays approximately linearly in a region between ∼15 and 25 mm from the magnet, resulting in a relatively consistent magnitude of force over this range. This region also coincides with a transition in the direction of force, with the magnet pulling  in the positive x-direction when z < 15 mm (as opposed to the negative x-direction, commensurate with the nominal direction of force, past this transition point). This spatial range is of interest as it approximately corresponds to the depth of the cerebral cortex, which contains vessel diameters typically between 2.5 and 40 micron [55,56], with flow velocities in the order of ∼0.5-1 mm s −1 [48,49]. Considering a normalized force of ∼10 T m −1 , and extrapolating from the simulations in section 3.4, an extremely high capture efficiency would be expected.
The field and force maps exhibited in figure 22 show that the most intense fields occur near to regions where the magnetization changes between segments in the array. Notably, the nominal direction of force shown in figure 22(d) points from the position of interest, towards the interface between the segments magnetized in the z-direction and the positive xdirection (gray and green segments respectively in figure 20). The change in direction in the force vector seen along the zaxis in figure 21(b) may then be understood as particles getting close enough to the magnet that they are more attracted to the local maximum in field adjacent to the interface between segments magnetized in the z-and negative x-directions (colour-coded gray and blue).

Conclusion
An optimization routine developed previously for the design of Halbach arrays consisting of cubic elements has been modified and expanded to generate arbitrarily shaped magnet arrays optimized to deliver magnetic force depending on a range of different design parameters. We have presented designs of optimized uniform magnet geometries and Halbach arrays, demonstrating how the performance of different arrangements varies as a function of the design parameters. The magn etic force applied by the arrays increases logarithmically with magnet volume, while the force emitted at the position of interest decreases almost exponentially as the position of interest gets further from the magnet. The number of allowed magnetization vectors is considered as a design parameter, and while using a greater variety of different magnetizations results in increased force output, it also leads to arrays that are more difficult to assemble, owing to repulsive dipole forces between neighbouring elements. A method to overcome this problem is considered, with simulations suggesting that introducing a small gap between repulsive segments can make optimized Halbach arrays significantly easier to assemble without causing a large compromise to the applied field and force. Simulations of magnetic microbeads under the influence of fluid flow in a 3 mm wide channel and a subset of optim ized array designs performed using COMSOL software suggest that a useful proportion of particles could be captured and retained at short range (5 mm) in mean fluid velocities of 100 mm s −1 or at further depths of 50 mm, when the velocity was 2.5 mm s −1 , depending on the choice of magnet. Finally, a design for a helmet magnet to apply an optimal magnetic force 50 mm deep inside the brain was generated to show the versatility of the optimization routine to address specific applications. For this design, particles tend to be most strongly attracted towards regions of the magnet where an interface exists between magnet segments. Based on the flow regime in the cerebral cortex, we suggest a high proportion of trapping in this region is feasible. Our examples of optim ized arrays show that, using the present optimization routine, magnet arrays can be designed and assessed for specific magn etic drug targeting applications once the required depth of targeting, direction of magnetic force, volume of magnet and the physiological features and flow regimes around the target have been accounted for.