Optimization and inverse modelling for MRI shielding applications

By Chris Riley

In this blog post I take a look at MRI shielding, both looking at the design of the shield and the inverse issue when there is a need to develop a suitable model for the magnet.

It’s true to say the first work in this area began in the 1950s but it was in the 1970s that MRI really took off. The two people who made this technology become what it is now, a common-place tool in hospitals – Paul Lauterbur of Stony Brook University in New York, and Peter Mansfield of University of Nottingham in the UK – were awarded the Nobel Prize for their work.

As far as electromagnetic modelling is concerned, we’re interested in the DC and RF coils and AC or periodic functioning coils inside the magnet, and these need to be modelled. In particular there is a requirement for a very highly homogeneous and very strong DC magnetic field over the image volume. The DC field aligns the spins of the atoms, and when an RF pulse is put into image volume area in the magnet, it excites the energy state of those atoms. The MRI system looks at the way the energy state relaxes back to its normal state, and you can tell what type of atoms are actually being dealt with. It’s localised by the system of gradient coils, which are pulsed coils, so that only one small part of the volume, at any particular time, is at exactly the right flux density to give the resonance when the RF pulse comes in.

An important thing to note is the higher the DC field that is used, the better the resolution that you actually get in the magnet. To achieve those high DC fields, certainly for a whole body scanner, the only viable technology to produce the DC field over the type of volume we’re talking about (e.g. maybe a spherical volume of half a metre in diameter, with very high homogeneity) is a set of superconducting coils. If you tried to do this with copper coils the resistive losses would be absolutely huge.

Here is a representative set of MRI coils:

1

This is the type of arrangement you see very commonly inside MRIs. They are all superconducting solenoid coils with many turns, and this arrangement of coils is trying to produce a homogenous volume in the middle. There are main DC coils which are making an axial field in one particular direction, and then there are also shielding coils around it, which are actually pushing flux in the opposite direction, trying to negate the stray field from the main DC coils.

We can see that even with those active shield coils we still have a field of around 7.5 Gauss, around two and half metres away from the magnet. To comply with legal requirements, that requires shielding. We need to do something to stop that level of field from being seen from that distance away from the magnet, usually because there are people present.

There is a legislation that is uniformly adopted throughout the world which says where humans are subject to a DC field, that field must be below 5 Gauss. That’s the first purpose of the shield, to try to divert the flux so it doesn’t leak out into the areas to where humans are present. The shield also serves a second purpose – it needs to shield the magnet itself from exterior magnetic sources. These are things like power cables. The reason it needs to be shielded is because we don’t want any stray fields from those cables or other magnetic sources getting into the centre of the magnet and perturbing the field. The shield has to serve in two directions. This even comes down to being able to shield from the presence of traffic. There are lots of shielding calculations done where it’s necessary to have the correct operation of the magnet whether there is a large metallic truck parked on the other side of the wall or not. But, importantly, we also do not want that shield to prevent the correct operation of the MRI magnet.

MRI magnets have a very homogenous central field, so a small perturbation in the return path of the field will change the homogeneity of the magnet. What’s important to the magnet manufacturers is that it doesn’t change by too much so that they can still be able to shim the magnet to produce the level of homogeneity that they require.

Of course, all these requirements also need to be done at a minimal cost, so basically with the minimum amount of material used in the shield. It may be easy to design a shield that works, but it may be very expensive.

How do we simulate this in Opera-3d? It’s an application that’s done with the non-linear magnetostatic analysis, which solves for the magnetic scalar potentials. It’s a very efficient calculation because it’s only solving for one variable at each node. We’ve introduced a thin sheet boundary condition, which is a 2d representation of the shield. Instead of actually trying to make a volume of the plates that are used to make up the shield materials, we can actually just represent this as a 2 dimensional surface.

The user specifies what material it is, and how thick that shield is, although there is no finite volume to the boundary condition. If we wanted to do a layered shield, made up of several plates, we could also give a packing factor. This is not 100% solid steel, this is perhaps 95% solid steel, because there are small gaps between the plates.  We’ve also introduced post-processing tools which are particularly useful for MRI designers.

Let’s take a look at a thin plate boundary condition first:

2

This is just the 2 dimensional surface. The user has said that on these surfaces surrounding the magnet a boundary condition called “Thin plate” will be placed. The only data that is actually necessary is the material name, which is steel. This will equate to a magnetic characteristic when we define the material properties. For the plate thickness in this particular case we’ve used a parametric value, which is just representing a particular number. Parameterisation is an important part of optimisation.

That’s the thin plate boundary condition. It certainly makes model building easier, and it usually makes analysis a lot faster as well.

In the post-processing we’ve got the ability to plot iso-contours of flux density:

3

We can see that this shield is not adequate at the moment. Most of the contour for the 5 Gauss line is outside the shield still. The post-processor allows us to look at this in more detail. We can see where the 5 Gauss contour finishes. We can look in a lot more detail at what the stray field is at any point in space.

I also mentioned that assessing the effect of the field on the homogeneity is important, and the standard way that MRI designers assess field homogeneity in their magnets is to look at the associated Legendre polynomial coefficients on a spherical surface, because that harmonic series is a solution to Laplace’s equation. Laplace’s equation is valid in a source free volume. There’s nothing, there’s no sources inside that sphere so we can find out the field anywhere inside that volume. Just looking at the harmonic coefficients is a good way of being able to compare designs, and whether we’ve improved the homogeneity, or reduced it by what we’ve actually done. Opera has a tool to be able to do this fit of the associated Legendre polynomial coefficients on the spherical surface.

Users are likely to need to optimise the design, and Opera includes an Optimizer module. This works with not just the magnetostatics simulation, but with any other type of simulation that you want to run with the Opera-2d or the Opera-3d software. It can be constrained or unconstrained. You can use any variable that you want to inside a model, as long as we can define a model in terms of variables. This could be the geometry we’re defining in terms of variables, e.g. for the thin plate boundary condition the thickness had been parameterised, so we could use that parameter as one of the variables in the optimisation. But we could also be using other things like material properties or excitation – anything that we can consider as a continuous variable.

The Opera-3d Modeller, which is the tool that we use to create models, allows us to perform parameterisation. This can also be done on CAD models – these can be imported into Opera from another package, and the Modeller can be used to make modifications to that model that will allow it to be parameterised. It’s multi-variable, and it’s also multi-objective, so we don’t just have to try to achieve one objective – we can achieve several objectives at the same time. So we could look at whether the shield was functioning properly, but also look at whether it is perturbing the fields inside the magnet.

It uses the Pareto optimisation technique to determine a family of best designs. Because we’re now dealing with multiple objectives, it will output more than one design, some which will be very good for one particular objective, some are very good in another objective. Hopefully there are some in the middle that are good on all objectives and typically that’s the type of design that people will want to use. Those objectives can be anything that the user can possibly evaluate within a post-processor. So the objectives are calculated using a script inside our post-processor, which enables the user to evaluate objectives. That can be as complex as you’d like it to be, it can even go and run Python code to calculate other objectives if the user wants to.

Here’s a simple example:

4

On the right side we’ve put a single plate. We’re currently looking at the contours of the flux density within that plate, which we want to try and optimise the shape of, to try and achieve certain objectives.

The first objective we want to achieve is to not use any more steel than we have to. The second objective is that we want to minimize the flux outside the shield, and the third objective is that we want to minimize perturbation of the field homogeneity. We can also add constraints to that: we also want to add that everywhere 10cm outside the shield is less than 5 Gauss, and we don’t want to see the central field of the magnet change by more than 15 Gauss, because we want to be able to shim it later.

The plate is 2m from the magnet centre and is defined by 3 dimensions – its length, its height and its thickness, and of course that thickness is applied to the 2d thin plate boundary condition. The user can specify in the Optimizer the range over which those variables will be searched. We’re defining the volume of variable space that we want to examine those values in. There can be constraints applied from that, by limiting the range. If we know that we can never accommodate a plate over 600cm high, then we can just limit the range we want to look at. Even though we may have used a lot of parameters inside the actual model, we don’t have to use all of them.

The objectives are evaluated using a script for the post-processor. It’s running two particular commands. It’s examining a field on a rectangular surface just outside the plate, 10 cm outside. It’s performing a map of the flux density, although that is not strictly necessary. That is just from a visual point of view. It’s also running the Legendre polynomial coefficient calculation on the axial component of the flux density on a sphere with a 25cm radius.

5

From those two commands we’re using results have gone into the system variables in the post-processor which we can then use to create our output values that will return to the Optimizer. I’m setting up a user constant called #intB which is the integral on that surface of component B. I’m also setting up another one to find out its maximum value. I’m also setting up export parameters for the harmonic coefficients as relative values to the central field, because that central field will change every time we change the shield design. We can run any commands like these to calculate complex expressions within the script.

I also want to set up another variable which gives me the difference between the central field and the value it’s currently set at. And finally I’ll add a quick calculation for the volume of the shield.

I don’t have to use everything so in this particular optimisation and I’ve only used 3 of the 6 quantities that are being calculated by the post-processor through that script as objectives. The Optimizer offers the opportunity to set objectives to minimize values, and also to maximise values.

6

This is what the Optimizer will produce. It looks at a series of models where it’s trying to produce the optimised designs. Because it’s multi-objective, it comes up with a family of designs. In this instance there are about 10 or 11 designs there in blue in which it satisfies the constraints and which are the best results. It is ranking the results received. Further down the list there will be some designs which have not satisfied the constraints.

We can see this when we look at the graphical presentation of these results…