Modeling Magnetic Particle Imaging for Dynamic Tracer Distributions

Magnetic Particle Imaging (MPI) is a promising tracer-based, functional medical imaging technique which measures the non-linear magnetization response of magnetic nanoparticles to a dynamic magnetic field. For image reconstruction, system matrices from time-consuming calibration scans are used predominantly. Finding modeled forward operators for magnetic particle imaging, which are able to compete with measured matrices in practice, is an ongoing topic of research. The existing models for magnetic particle imaging are by design not suitable for arbitrary dynamic tracer concentrations. Neither modeled nor measured system matrices account for changes in the concentration during a single scanning cycle. In this paper we present a new MPI forward model for dynamic concentrations. A static model will be introduced briefly, followed by the changes due to the dynamic behavior of the tracer concentration. Furthermore, the relevance of this new extended model is examined by investigating the influence of the extension and example reconstructions with the new and the standard model.


Introduction
Magnetic Particle Imaging (MPI) is a relatively new medical imaging modality invented by Weizenecker and Gleich in 2005 [9]. In this tomographic imaging technique, the non-linear magnetization response of the superparamagnetic tracer material to an external magnetic field induces a potential in the receive coils of the scanner. The spatial distribution of the magnetic particles is reconstructed from these measurements. MPI allows for a rapid data acquisition with high temporal resolution which makes it a promising imaging device for different imaging applications, see [21] for an overview. In many of these applications, visualization of tracer dynamics is highly relevant, such as physiological diagnosis like stroke detection [24], visualization of blood flow [31] or localization of medical instruments in vascular interventions [14]. The MPI forward operator can be described by model-or measurement-based approaches [13]. In a measurement-based approach the forward operator is represented by a calibration scan [27,30]. Therefore, the signal generated by a delta sample of tracer material is measured for a finite number of spatial positions. The modeling approach describes the measurement process with physical laws [22]. Unfortunately, models usually idealize the physical setting to limit the complexity of the model. These simplifications can lead to large modeling errors and give reasons for the time consuming measurement approach being still dominant in practice. The state of the art model underlies the assumption of a (nearly) static concentration during the signal acquisition. This assumption is not always fulfilled. MPI is able to visualize the distribution of a liquid tracer. It can accumulate, dissipate or move e.g. with the blood flow. The behavior of the particle concentration is not static in these cases. Also time-series measurements imply a dynamic tracer distribution such that the static model is only true for piecewise constant concentrations. The same problem is valid for the measurement approach. Since the delta sample is static during each cycle of the calibration scan, the measured system matrix does not cover dynamic behavior. Currently, the only way to reconstruct non-periodic dynamic concentrations is to reconstruct a time-series of images under the assumption of static behavior during the scan [28,14]. Reconstruction of periodic dynamics in magnetic particle imaging is investigated in [8] in order to reduce of artifacts induced by cardiac-or respiratory motion in multi-patch MPI. The authors use the measurement-based approach and assume limits on the velocity and periodicity of the motion. By rearranging measurements from the same motion phase into virtual frames, dynamic tracer distributions can be reconstructed by static reconstructions from the virtual frames. The model-based approach gives rise to various directions of research covering all components of the signal generation chain and analyses of the models. One of these directions is modeling of the magnetic behavior of magnetic nanoparticles which was studied by Kluth [17,19] and Weizenecker [29]. The most frequently used magnetization model is the Langevin-or equilibrium model, which is also the basis for the derivations in the following articles. The equilibrium model does not respect magnetic relaxation effects. In [17], the model is extended for different kinds of relaxation. The author presents forward models incorporating either Brownian rotation or Néel relaxation in the cases of mono-and polydisperse tracers under the assumption of single domain particles with uniaxial anisotropy. Based on the equilibrium model the authors of [26] derive analytical reconstruction formulae as well as numerical reconstruction schemes for two-and three-dimensional MPI. They examine and compare the ill-posedness of the reconstruction problem for different dimensions. A mathematical analysis of the 1D model is provided by Erb, Weinmann et al. [7]. They investigate properties like the ill-posedness and discover an exponential singular value decay of the reconstruction problem. Goodwill and Conolly [10] follow the X-space approach. They consider the dependence of the spatial position of the field free point (FFP), which is the time-dependent volume of vanishing magnetic field strength, and the drive field of the scanner. As a result the forward operator of 1D MPI can be identified as a convolution operator. The authors extend their approach to multiple dimensions in [11]. The more analytically focused article by Maass and Mertins deals with closed-form expressions for the Fourier transform of the system function for multiple dimensions [25]. The system function is related to tensor products of Chebychev polynomials of the second kind and tensor products of Bessel functions. This result might allow for analytical insights into the system function and more efficient reconstruction techniques in the future.
Another common model simplification is the assumption of ideal magnetic fields. In practice magnetic fields can be distorted which is influencing spatial signal encoding. In [5], the authors use spherical harmonics to achieve more realistic representations of magnetic fields for FFP and field free line (FFL) scanners. A 3D forward model for non-ideal magnetic fields which can be reconstructed with the algebraic reconstruction technique (ART) is presented. Artifacts caused by distorted magnetic fields have also been investigated in the context of the X-space approach [32]. These distortions especially affect multi-patch MPI, since distortions increase with the distance from the center. In the case of measurementbased reconstruction, compensation methods for displacement artifacts in multi-patch scans are studied in [4,3]. In this paper we present an extended MPI forward model for dynamic tracer distributions, in the discrete and continuous case, both in time-and Fourier domain. While the initial theoretical setup presented in Sec.3 is identical to the one in [8], our model is not limited to periodic motion and covers dynamic tracer distributions with high velocities. Furthermore, we provide simulation experiments to show the importance of the extension relative to the tracer dynamics and the impact on reconstruction quality compared to the static model. The presented approach is of special interest for blood flow measurements [16] because the speed of the motion is part of the model and can be reconstructed simultaneously. The remainder of this paper starts with a brief introduction to the principles of an ideal MPI system in Sec.2 and is followed by presenting the standard modeling approach in Sec.2.1 which we will extend to arbitrary dynamic tracer distributions in Sec.3. Based on the Langevin model, FFP scanners and Lissajous trajectories we investigate the influence of the extension to the signal for different kinds of dynamics in Sec.4, while we compare reconstructions of simulated dynamic measurements with the standard and the extended model in Sec.5. We close with a discussion of the results in Sec.6.
2 Basic principles of magnetic particle imaging The aim of magnetic particle imaging is the reconstruction of the multi-dimensional spatial concentration of the particles. Spatial encoding of the signal is realized by applying a spatially and temporally varying magnetic field H ∈ L 2 (Ω × R + , R 3 ) where Ω ⊂ R 3 denotes the field of view (FOV). The magnetic field consists of a spatially inhomogeneous selection field H S and a temporally varying drive field H D . The selection field is a gradient field that has a point of zero field strength in the center, the so-called field free point (FFP), and a linearly increasing field strength to the periphery (see Fig.1a). The drive field H D (t) = [a l sin(2π f l t + ϕ l )] l=1,...,3 is spatially constant but changes its magnetization over time according to a sine function in each dimension. It has three parameters per dimension, the amplitude a l ∈ R determining the size of the field of view, the frequency f l ∈ R defining the density of the scan trajectory and a phase shift ϕ l ∈ R setting the starting point of the scan trajectory. Choosing the parameters appropriately, the overlay of H S and H D forms the field H which has a FFP moving through the volume of interest along a so-called Lissajous curve (see Fig.1b). Superparamagnetic means that the magnetic nanoparticles behave like tiny magnets, while an external magnetic field is applied. They have their own magnetic moments which are larger than their atomic moments. There is no remanent magnetization after the applied magnetic field is removed [2].  The magnetic moment of the particles responds to temporal changes of magnetic fields. There are different models describing the magnetic behavior of the particles which where studied in [17,19,29]. In Sec.4 and 5 the Langevin or equilibrium model is used but could be replaced by more complex models. When the field free point moves over a position r, it causes a change in the mean magnetic moment at this location. The magnitude of the magnetization is proportional to the tracer concentration c : Ω → R + and the mean magnetic moment m : Ω × R + → R 3 . The Sobolev space denotes the space of L 2 -functions whose first weak derivatives are also functions in L 2 on the domain D ⊂ R d with γ ∈ N d being a multi-index with |γ| = ∑ d i=1 γ i , see [1]. The change in magnetization induces a current in the receive coils of the scanner. Due to the construction of the magnetic field, a measured voltage at time point t can be connected to a magnetization change and thus to a certain concentration at a position r. The measurement process is described by a forward model in the following section.

MPI forward model
The static forward model describes the magetic particle imaging process in time domain [20]. The sensitivity of the receive coils p is multiplied with the permeability constant µ 0 and the change in magnetization which is caused by the superparamagnetic particles and the applied magnetic field H from the send-coils.
Assuming that the signal generated by the excitation field is removed by a filter yields As defined in the previous section, the excitation field is multi-dimensional which means that the system function S(r,t) maps to R 3 , thus u(t) ∈ R 3 is a voltage vector, where each value is measured by the respective receive coil. In the discrete models in the remainder of this paper we will refer to a single component of u(t) since the computations are analogous for all channels.
Discretization The formulation of discrete forward models is motivated by measured system matrices and the use of numerical reconstruction methods. Therefore, we use a basis An intuitive choice are piecewise constant basis functions on equisized, pairwise disjoint quadratic or cubic domains as they are a reasonable representation of both the pixels or voxels in an image and the delta probe used for the calibration scans.
Using the basis functions, we obtain piecewise constant approximations of the concentration and system functioñ and T c being the repetition time for one Lissajous cycle. Instertingc in (3) yields the following discrete forward problem It can also be written as a matrix vector multiplication of a concentration vector and the system matrix S Reconstructing the concentration vector c from a given measurement vector u is a classic inverse problem. In [18,26], it was shown that the multidimensional MPI reconstruction problem is severely ill-posed. Thus, computing a stable and unique solution requires regularization. Two common regularization methods in MPI are Tikhonov-and iterative regularization. The former defines a Tikhonov functional by adding a penalty term with a regularization parameter. The resulting minimization problem can then be solved by the Kaczmarz algorithm or other iterative schemes adapted to the applied regularization term. The latter option regularizes the iterative method directly by choosing a maximum number of iterations. In both cases, the standard iterative method used in MPI is the Kaczmarz algorithm [15]. The information from reconstructions of several channels can be combined to improve image quality.

Dynamic forward model
The forward model presented in the preceding section underlies the assumption of a (nearly) static concentration during the signal acquisition which might be violated in case of dynamically changing tracer distributions. MPI is able to visualize the distribution of a liquid tracer which can accumulate, dissipate or move e.g. with the blood flow. In these situations, the behavior of the particle concentration is clearly not static.
In practice oftentimes measured system matrices are used for MPI reconstruction. These matrices are the results of calibration scans which measure the induced voltage of a delta sample during a scanning cycle for each spatial position. This approach yields good results for static concentrations because the matrices also incorporate the transfer function of the system. Since the delta sample is static during the complete cycle the measured system matrix does not cover dynamic behavior.
In order to adapt the model to dynamic tracer concentrations, the magnetization function (1) is modified such that it contains a time-dependent concentration Thus, the forward model (2) changes to Assuming a constant coil sensitivity p and that the signal generated by the excitation field is removed by a filter results in the dynamic forward model It describes a measurement u : R + → R 3 in time domain and contains a sum of two system functions S 1 , S 2 : Ω × R + → R 3 multiplied with the tracer concentration c : Ω × R + → R + and its time derivative. The derivatives ∂ c ∂ t , S 1 = ∂m ∂ t and the measurement u are L 2 -functions because the concentration c and S 2 =m are in the Sobolev space H 1 (Ω × R + ).
Dynamic forward model in frequency domain MPI measurements are usually given in frequency domain. Due to the time dependence of the concentration the static model The measurement in frequency spaceû : R + → C 3 and the derivatives ∂ c ∂ t : Ω × R + → C, The convolution is only applied to the frequency components.
Discretization Using the same pixel-basis {φ i } i=1,...,R ⊂ L 2 for discretization as in Sec.2.1 yields the following representation of a piecewise constant dynamic concentratioñ and analogously for the derivative ∂ c ∂ t and system functions S 1 and S 2 . Together with the time sampling points {t j } j=1,...,n T from Sec.2.1, we get the discretized dynamic forward problem For measurements with F ≥ 1 cycles the time sampling for the measurement and concentration changes to {τ j } j=1,...,Fn T with τ j = ( j − 1)FT c /(Fn T − 1) while the system functions are evaluated at t j mod n T . Eq. (12) is no longer a matrix-vector multiplication as in (5) but a sum of element-wise matrix multiplications ..,n , with matrices A, B ∈ R n×m . In frequency space the same approach yields the following discrete forward problem Again the convolution is only applied to the frequency components, i.e. the respective matrix columns. Note that a frequency domain reconstruction computesĉ. To see the behavior of the concentration in time, the inverse Fourier transform needs to be applied. Reconstruction becomes a deconvolution problem in frequency space. A typical solution approach for this ill-posed inverse problem is to make use of the convolution theorem of the Fourier transform which in this case results in time domain reconstruction. The dynamic model (8) is also mentioned in [8] but followed by strong restrictions of the dynamics such that there are no further consequences in the reconstruction process. In contrast, the models proposed in this section are valid for a broad range of dynamics, e.g. rapid changes or non-periodic behavior. The tracer distribution is required to be differentiable in time and integrable in space.

Relevance of the dynamic model
As mentioned in Sec.2.1, the concentration is usually assumed to be constant. The time derivative of the concentration would therefore be nearly zero and the second summand of the extended model (9) would thus be small such that it can be neglected. We investigate the structure of S 2 in comparison to S 1 and consider a set of simulated dynamic concentrations to survey whether neglecting the second term in the new model may Langevin model In the Langevin model the particles are assumed to be in thermal equilibrium and the applied magnetic fields to be static. The mean magnetic moment at spatial position r and time point t is given bȳ with L α,β : R → R being the Langevin function with α, β ∈ R being particle dependent parameters.
Eq.(11) shows a sum of two convolutions. In a first step, we are interested in the shape of the convolution kernels. Therefore, we compute which are shown in Fig.2 together with an approximation of their convex hulls. The approximation of the convex hull was calculated by determining and connecting the maximum values within the next 15 frequency steps to include all peaks of the function. Both matrices exhibit a similar structure. The peaks have the same distances (≈ 15 frequency steps) and the convex hull (orange line) has a similar shape. The full-width-at-half-maximum (FWHM) of the convex hull is the same (≈ 33 frequency steps) while the maximum of the second system matrixŜ 2 is 10 4 smaller than the maximum of the first system matrixŜ 1 . Thus, on a first view, the assumption that the second term S 2 in the dynamic model is negligible might be reasonable. In a second step, we looked at four types of dynamic concentrations during one cycle for a single voxel. The plots in Fig.3 show the concentration over the scan time, its time derivative and the respective Fourier transforms for each example concentration. Example concentration 1 is depicted in Fig.3a which shows one peak at the beginning of the scan. The tracer is flowing through the voxel for a short period of time. This could be a small tracer bolus moving fast through the volume of the voxel. In the second example, shown in Fig.3b, the concentration increases strongly in the beginning, remains constant for a short period of time and decreases again. The tracer flows through the voxel for a longer period of time. Example 2 represents a larger bolus moving fast through the voxel. Example 3, shown in Fig.3c, shows a slow increase and decrease of the concentration. This represents a slowly moving small bolus. Example 4 is a periodic version of the first example. Fig.3d shows two peaks within the scan time. The tracer flows two times through the voxel with a high velocity. This represents a small bolus with fast periodic motion.
Looking at the Fourier transformations shows that the maximal absolute values of the Fourier transformed concentrationsĉ are about 10 4 smaller than the maximal absolute values of the Fourier transformed time derivatives dc dt for all 4 examples. This demonstrates that for these dynamic concentrations the magnitude of the two summands of the new dynamic model is the same. The imaging process in Fourier space is a convolution of the Fourier transformed system matrices with the Fourier transformed concentration and its time derivative. Thus, the concentration is smoothed by the system matrix. The kernels S 1 and S 2 have the same width meaning that the concentration and its derivative are smoothed equally. To further examine these effects, we split up the discrete forward model such that the signal is now the sum of A and B, where A denotes the signal component generated by the first system matrix S 1 and B the signal component generated by the second system matrix S 2 .
The convolution of the frequency components of system matrix 1 and the tracer distribution is named a and the convolution of system matrix 2 with the derivative of the concentration is named b. Fig. 4: Absolute values of the convolution of the system matrices with the concentration and its derivative. Each curve shows the frequency amplitudes for one of the 19 2 voxels.
Using the dynamic example concentration 3 shown in Fig.3c, Fig.4 shows a and b for each voxel. As expected, one can see that the shape and the maximum values of both terms are similar. Both plots show maximum values of about 1.2 · 10 −4 . Fig.5 shows the plots of A and B for the frequencies k ∈ [0.08, 1.25]MHz. Again, both terms have the same order of magnitude. For A the frequencies with high amplitudes have a small variance while for B high amplitudes can be observed in the whole frequency range. This shows that even for this example with slower dynamics the second component of the forward model has a significant impact on the signal. Thus, the second summand in Eq.(9) should not be generally neglected for dynamic concentrations.

Comparing reconstructions with the dynamic and static model
The challenge in solving the dynamic inverse problem (12) is the high number of degrees of freedom. We therefore use a minimalist setup with a grid of 3 × 3 × 1 voxels. We use two computational phantoms to simulate an MPI measurement with the dynamic model and reconstruct it with the dynamic and the static model. They are named one-peak phantom and (a) Signal part A = ∑ R i=1Ŝ1 (r i ,k j ) * ĉ(r i ,k j ) generated by the first system matrix three-peak phantom and their spatial setup can be seen in Fig.6.  Fig. 6: Spatial setup of the One-and three-peak phantom. They consist of 3 × 3 × 1 voxels indexed from 1 to 9.

Parameterized concentration curve
To reduce the degrees of freedom and computational cost, the tracer concentration of the one-peak phantom is described by parametric curves ..,R , M ∈ N and basis functions ψ m which are cubic B-splines. This means that for each voxel r i there is a set of coefficients b m which together with the spline basis form a continuous concentration curve in time. Consequently, we assume that the concentration is twice differentiable with respect to t which is a stronger condition than previously assumed in the dynamic model. Cubic B-spline curves are well suited to model the dynamics of the magnetic tracer. In [12], the authors deal with the reconstruction of spatiotemporal tracer distributions in Single Photon Emission Computed Tomography (SPECT). Cubic B-spline curves are used to describe and reconstruct the dynamic distribution of the radioactive tracer from gated cardiac SPECT sequences. We generate three variants of the one-peak phantom 1F, 2F and 4F. For one-peak phantom 1F the coefficients for all R = 9 voxels except r 5 are zero. The concentration is non-zero within the scan time of one frame. Fig.7a shows the development of the tracer distribution for the total scan time where each curve describes the concentration within one voxel. The plot shows a peak at t = 0.4128ms with a concentration of 2.67 for voxel r 5 . Versions 2F and 4F differ only in the width of the concentration peak. The concentration peak for r 5 lasts for the scan time of 2 frames in version 2F and 4 frames in version 4F (see Fig.7b). The three variants can be related to boluses with different velocities. The bolus in version 1F is twice as fast as in 2F and four times faster than in version 4F. Measurements with 4 frames which are each sampled at 408 time points and the dynamic forward model with S 1 = ∂m ∂ t and S 2 =m are simulated according to (12). The parameters used for the simulations are listed in Tab.1 and Tab.2. There is no transition time in between the frames. We reconstruct the concentrations by minimization with respect to the parameter set λ B so that we get continuous concentration curves. The solution set is restricted to parametric (a) Reconstruction of 1F with the dynamic model using S 1 and S 2 .
(b) Reconstruction of 1F using only S 1 .
(c) Reconstruction of 2F with the dynamic model using S 1 and S 2 .
(d) Reconstruction of 2F using only S 1 .
(e) Reconstruction of 4F with the dynamic model using S 1 and S 2 .
(f) Reconstruction of 4F using only S 1 . Fig. 8: Measurements of the dynamic one-peak phantoms 1F, 2F and 4F are simulated with the dynamic forward model (12). They are reconstructed with either both S 1 and S 2 (left) or only S 1 (right). All plots show averages of x-and y-channel reconstructions. The dashed lines outline the true concentration in voxel r 5 and the vertical grid lines mark the start and end of frames. spline curves in L 2 (R 3 ) × C 2 (R) which is an implicit regularization. The curves are reconstructed with two different settings. In the first experiment both matrices are used for reconstruction which corresponds to minimizing with u ∈ R n T . The problem is minimized with 200 iterations of a conjugate gradient algorithm and no further regularization. Fig.8a shows the average of reconstructions of the x-and y-channel of one-peak phantom 1F. The peak for voxel r 5 is located at t = 0.4448ms with a concentration of 2.96 which is very close to ground truth. In the same period of time also the concentration of the remaining voxels is non-zero. The peaks of the voxels with even indices have concentration values of about 0.9 and the peaks of the voxels with odd indices have even smaller values of about 0.4. Even if these voxels have a non-zero concentration, it is significantly lower than the value of r 5 , so that we can expect sufficient contrast in the reconstructed images. The values for the off-diagonal voxels (voxels with even indices, cf. Fig.6) show higher concentration values than the ones on the diagonal. The x-channel reconstruction locates the concentration correctly in x and the y-channel reconstruction locates the concentration correctly in y. Thus, the off-diagonal voxels are masked by the high concentration in the central voxel.
In the next experiment the same measurement is reconstructed only with S 1 which corresponds to minimizing with u ∈ R n T . The problem is again minimized with 200 iterations of a conjugate gradient algorithm and no further regularization. The result for the average of x-and y-channel reconstructions is shown in Fig.8b. The reconstruction shows a peak for voxel r 5 at t = 0.3808ms which is close to the ground truth but with a significantly smaller concentration of 0.85. Again there are concentration peaks for all remaining voxels with values of about 0.4. This means that the reconstructed images will show reduced contrast. And the true concentration is underestimated. To get a more intuitive impression of the impact of the discussed curves on the reconstruction quality, Fig.9 shows a frame of the phantom and the two reconstructed time-series at the time point of the maximum concentration (t = 0.4128ms). Looking at the image of the first experiment in Fig.9b one can observe a good contrast and slightly higher concentration values for the off-diagonal voxels. As discussed above, the image from the second experiment shown in Fig.9c exhibits poor contrast and significantly lower concentration values compared to the phantom. Version 2F and 4F of the one-peak phantom are reconstructed analogously to version 1F. While also for version 2F reconstructed with both matrices the concentration peaks show the correct location and 83% of the true amplitude (see Fig.8c), the peak in the reconstruction using only S 1 is less than 50% of ground truth (see Fig.8d). The reconstructions for version 4F with and without S 2 are shown in Fig.8e and 8f. The location is correct for both reconstruction methods and the amplitude of the peak for r 5 reaches 87% of the ground truth for reconstruction with both matrices and 73% for using only S 1 . Also the concentrations for the remaining voxels are sufficiently low in both cases. Looking at a frame of one-peak phantom 4F and the two reconstructed time-series at the time point of the maximum concentration (t = 1.304ms) in Fig.10, one can see that the contrast in the reconstruction without S 2 (see. Fig.10c) is improved compared to version 1F and almost comparable to the reconstruction using S 1 and S 2 (see. Fig.10b). In order to get an impression on the strength of the dynamics in one-peak phantom 4F, we relate these values to a 2×2×2 mm 3 bolus with constant concentration c max moved through a 2 × 2 × 2 mm 3 voxel with a constant velocity v (see Fig.11). The one-peak phantom has a maximum concentration of c max = 2.667 and a maximum time derivative ofċ max = 3065. The average change rate yields a bolus velocity of v av = 2 · 10 −3 /(2 · T c ) = 1.53 m/s and the maximum change rate results in a velocity of v max = 2 · 10 −3 /(c max /ċ max ) = 2.3 m/s. Thus, we state v dyn = 1.53 m/s as a preliminary threshold velocity. For reconstructions with an average flow v > v dyn the dynamic model will improve the reconstruction quality in comparison to the static model.

Frame-by-frame reconstruction
Another way to investigate the impact of the new model is to reconstruct a dynamic measurement frame-by-frame with the assumption of a static concentration within each frame. Therefore we use the three-peak phantom. As the one-peak phantom, the tracer concentration of the three-peak phantom is described by parametric curves c(r i ,t) = ∑ m b m,r i ψ m (t) = c(λ B (r i ),t) with parameter set λ B and basis functions ψ m which are cubic B-splines meaning that for each voxel r i there is a set of coefficients b m which together with the spline basis form a continuous concentration curve in time.
For the three-peak phantom only the coefficients for the voxels r 4 , r 5 and r 6 are non-zero. The tracer distribution during the total scan time is shown in Fig.12a where each curve describes the concentration within one voxel. There is a concentration peak of 6.67 for voxel r 4 , r 5 and r 6 . The peaks are shifted in time, such that this dynamic can be seen as an object or tracer bolus that moves from voxel r 4 to voxel r 6 considering the location of the voxels in Fig.6. The peaks are located in the scan time of frame 3, 4 and 5 and have a temporal width of about 4 frames. The concentration of the remaining voxels is zero. A measurement with F = 10 frames which are each sampled at 408 time points and the dynamic matrix model (12) with S 1 and S 2 is simulated. The parameters used for the simulation can be found in Tab.1 and Tab.2. The dynamic tracer distribution is reconstructed with two different settings. The first one uses information about the tracer dynamics from the reconstructions of previous frames and the second one reconstructs each frame independently. In fact the reconstructions are piecewise constant functions over time. For better comparison the results depicted in Fig.12 show linear interpolations of the static reconstructions of 10 frames. In the first experiment both matrices are used for reconstruction of each frame while the time is the divided difference of the concentration vector of the current and the preceding frame. This corresponds to minimizing Fig.12b shows the average of x-and y-channel reconstructions which were reconstructed in time domain with 100 iterations of a gradient descent algorithm and no further regularization. It can be seen that the peaks are correctly located in frame 3, 4 and 5. The amplitude of the peaks is slightly lower than the ground truth and decreasing, 5.41 for r 4 , 4.96 for r 5 and 4.82 for r 6 . There is a non-zero concentration for the remaining voxels in the first 5 frames of less than 0.5. So the reconstructed images will exhibit sufficient contrast.
In the next experiment the same measurement is reconstructed using only S 1 , i.e. minimizing with u f ∈ R n T , c f ∈ R R , S 1 ∈ R n T ×R in time domain with 100 iterations of a gradient descent algorithm and no further regularization. The result is shown in Fig.12c. Again the peaks are located correctly in frame 3, 4 and 5. The amplitudes 5.47 for r 4 , 5.07 for r 5 and 5.49 for r 6 are also slightly lower than in the phantom and differ less than in the first experiment. The remaining voxels show non-zero concentrations up to 1.2 being more than twice as high as for the first experiment.
(a) The concentration of the three-peak phantom changes in time only in voxel r 4 , r 5 and r 6 . In the remaining voxels the concentration is zero. The time-shifted concentration peaks form a motion from r 4 to r 6 .
(b) Average of x-and y-channel frameby-frame reconstructions with the dynamic model using S 1 , S 2 (c) Average of x-and y-channel frame-byframe reconstructions with the static model using only S 1 Fig. 12: A measurement of the dynamic three-peak phantom is simulated with the dynamic forward model (12). Each frame is reconstructed separately assuming a static tracer distribution within each frame. The frames are reconstructed with the dynamic and the static model.

Discussion and conclusion
We introduced a new extended forward model for dynamic magnetic particle imaging. It was shown that the standard forward model does not account for dynamic tracer distributions which is corrected by the extended model presented in this paper.
One of the main differences is that the dynamic model contains a second summand with a second system matrix. For different kinds of dynamic concentrations the two summands have been examined. The order of magnitude of the summands is the same for the chosen dynamic examples. This emphasizes the importance of the new model for dynamic tracer distributions. Furthermore, we simulated measurements from dynamic concentrations with the extended model and reconstructed them with both the dynamic and the static model. In the experiments in Sec.5.1 three simple phantoms with different change rates are examined. For onepeak 4F, the phantom with the lowest change rates, the static approach using only one system matrix provided an acceptable reconstruction quality. For the phantoms with higher change rates, one-peak phantom 2F and 1F, the static approach resulted in reconstructions with low contrast and significantly lower amplitudes than ground truth while the dynamic approach performed well on all three phantoms. While a quantitative study of this is beyond the scope of this article, we can state that for higher change rates than in one-peak phantom 4F the dynamic model should be considered for reconstruction. The presented dynamic model is more general than existing approaches for dynamic concentration reconstruction as it is not limited to periodic motion and can be applied to motions with high velocities. While in this paper the equilibrium model is used, the dynamic model allows to incorporate more advanced magnetization models which could improve the reconstruction quality in the future. The reconstruction approach with parametric concentration curves features an implicit dynamic regularization. Additional spatial or temporal regularization, like sparsity in time and space, can be included easily. Moreover, the model allows for joint reconstructions of the particle concentration and its time-derivative which might be of special interest for blood-flow diagnostics. It remains future research to develop new reconstruction techniques for dynamic tracer distributions based on this model and extend it to multi-patch imaging sequences for larger volumes. Furthermore, the methods need to be evaluated for simulations with phantoms of realistic size and physical phantom measurements. Practice-oriented scenarios might require improved minimization schemes for reconstruction. The dynamic model might also be combined with measurement-based reconstruction. While using a calibration scan for S 1 , the second system matrix can be modeled and corrected with the transfer function. Alternatively, S 2 might be learned from its measured time-derivative S 1 . In addition to our preliminary statement, a quantitative study of phantoms with different velocities is required for a more precise proposition about when the dynamic model is necessary based on the level of dynamics and the desired reconstruction quality. A further theoretical research direction is the analysis of features like the ill-posedness of the dynamic reconstruction problem.

Declarations
Funding No funds, grants, or other support was received.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.

Availability of Data
The data sets generated during the experiments are available from the corresponding author on reasonable request.