Prof Piotr Cysewski

Ph.D Piotr Cysewski
email: piotr.cysewski@cm.umk.pl

Non-additive interactions of nucleobases in model dinucleotide steps occurring in B-DNA crystals

Keywords: many-body, intermolecular interactions, dinucleotide steps, 8-oxoguanine, B-DNA Many-body  Intermolecular interactions  Dinucleotide steps  8-oxoguanine

Non-additivity of base-base interactions in all ten possible model dinucleotide steps were analyzed on MP2/aug-cc-pvDZ quantum chemistry level. Selected conformations of four nucleobases exactly matched to ones occurring in B-DNA crystals. In most of 162 analyzed tetramers both threeand four-body contributions are negligible except of d(GpG) steps. However, in these dinucleotides both contributions are always of opposite signs and in all cases the sum of all non-additive part of intermolecular interactions do not exceed 2.6 kcal/mol. This stands for less than 5% of the overall binding energy of dinucleotide steps. Also replacements of guanine with 8-oxoguanine in d(GpG) systems introduces non-additivity of the same magnitude as for canonical dinucleotides. It is observed linear relationships between values of the total binding energy obtained in the tetramer basis set and estimated energy exclusively in dimers basis sets with assumption of pairwise additivities. For all analyzed dinucleotides steps there are also linear correlations between amount of non-additive contributions and magnitude of pairs interactions. Based on differences in electrostatic contribution to the total binding energy of four nucleobases and polarity of dinucleotide steps three distinct classes of dinucleotide steps were identified.

Response to Reviewers: 1) I recommend to do a careful proofreading, there are many typos and some sentences are difficult to understand. Thank you. Careful proofreading was applied and text was modified accordingly.

2) Page 3, Results and discussion Alltogether 162 dinucleotide steps were analyzed. The number of individual steps is variable, why is it so? Why there was 48 d(GpG) steps, but only 10 d(GpC) steps? What were the PDB codes of structures from which steps were chosen, and which steps (chain ID, res ID) were chosen? What was the resolution of these structures? This information could go to the Supplementary material.

In supplementary material details related to analyzed structures are provided now. The differences in dinucleotide populations have strictly technical reason. Since in my previous work accepted for publication in Int.J.Quantum Chem. (

overrepresented. Since it constitutes its own class I do not think that such overrepresentation changes papers conclusions. On the other hand it is not necessary extending the number of other dinucleotide steps since the non-additivity is quite small and they were selected from broad range of data according to IIE and structural parameters diversities.

3) Page 3, Results and discussion

The numbers of steps are large enough to perform some statistical analysis of resulting energies. E.g. ANOVA could be performed to test whether the mean energies in dinucleotide steps differ significantly each from other. Also, is the distribution of energies for individual steps Gaussian or not? If not, some non-parametric testing procedures could be used. In general, the number of data seems to me to be large enough to do their thorough statistical analysis, which would considerably improve the quality of the results, maybe even leading to some yet unnoticed conclusions.

In fact I did statistical analysis but intentionally data were not discussed in details. The differences in the median values between all pairs of analyzed dinucleotide classes are greater than would be expected by chance (P = <0.001). Hence all three dinucleotide distributions are statistically different. The normality test failed only for the first class comprising d(GpG) dinucleotides. The other two sets passed both normality (P = 0.134) and equal variance tests (P = 0.096). I decided not including statistical analysis since the main goal of this paper was the answer to the question whether additivity is to be taken as important contribution to dinucleotide steps energetics. Thus, representative dinucleotide steps were taken into consideration and classified according to some common properties. Broad diversity of energy asymmetry coefficient was chosen as the selection criterion of dinucleotide steps. This ensured covering of the whole range of possible structures with minimized cost of computations. The presented quantification of non-additivities opens the door for more precise and broad characteristics of all possible di-, tri-and oligo-nucleotides. Since I have completed much more computations and I have now much extended point of view it seems to me that statistical analysis would be irrelevant in the context of this paper. Anyway I greatly appreciate this comment of the Referee.

4) Page 3, Method section

The MP2/avdz preoptimized geometries of monomers were used. I wonder what is the influence to the nonadditivy term if hydrogens were reoptimized within the base-pair? It would be worth evaluating this for some selected steps. I am asking this as in my experience there exist surprisingly not that small influence of hydrogen positions on the total stacking energies.

Thank you for this suggestion. This is really interesting subject, but I am afraid it is outside of the scope of this paper. 5) Page 3, Method section >From the Methods section it is not clear if the Cs symmetry was imposed during the monomer optimization. If yes, what is the influence to the nonadditivy term of keeping NH2 group, which is known to be pyramidalized, in plane? I am aware that non-planarity of amino groups is a fact. However according to my experience this is quite small effect and I do not think it is crucial for anon-additivity. In fact full geometry optimization is required for proper evaluation of pyramidalization of NH2 exocyclic groups. Again although this aspect is interesting I must emphasize that it not the main focus of my interest but will surely be considered in forthcoming investigations.

6) Page 4

The point BD0015 (3-body term of +24.2 kcal/mol and 4-bod of -21.2 kcal/mol) is not present in the Figure 2. This is of course related to the applied scale. I have extended the sale for showing all points. Thank you for this prompt, of course I did it with pleasure. 10) I wonder whether the sign of the nonadditivity term doesn't depend on the level of theory used? Maybe if calculated at the very high level of theory (CCSD(T) complete basis set limit) instead MP2/avdz, the conclusions about the attractive/repulsive nature could be different. It i clear that such high level of theory is not doable, however some discussion would be appropriate. I am afraid that at this moment the CCSD(T) level is beyond my computations capabilities for dinucleotide steps. Hence I only added warnings in the final section of my paper. 11) Also, MP2/avdz does not include dispersion non-additivity, this should be discussed (with appropriate references) in the text. Thank you. Suitable comments were added.

Introduction

Many-particle systems as immanent nature of physical world stand for emergence of new properties that cannot be derived just as simple adding up features of constituting elements.

Non-negligible contributions of many-body terms to intermolecular interactions in molecular clusters are well known [1][2][3][4] and this aspect significantly complicates assessments of stabilization energies in multi-molecular systems. For example three-body terms have significant contribution in water clusters reaching even 30% to the interaction energy [2].

Fortunately, four-body and higher order terms are usually negligible [1][2][3][4][5][6][7][8]. This is of particular importance since almost all of molecular modeling procedures ignores non-additive contributions to the total energy of molecular systems. Typical example is the contemporary molecular dynamics simulations for biomolecules. Size of these kind of systems enable merely for application of force field parameterization inherently neglecting not only quantum effects but also all non-additivities of interacting systems. Implementing hybrid QM/MM methodology introduce improvement only to some extent. Furthermore, corrections for electron correlation are important contributions to amount of individual n-body terms frequently imposing increase of their magnitudes. On the other hand many-body contributions can substantially vary with the relative orientation of monomers [5,6] and can originate from distinct types of interactions. For example in case of water molecules energies of hydrogen bonding are two order stronger than for van der Waals interactions in noble gas clusters.

Although for electrically neutral and non-polar systems non-additivity is typically small [7] for polar systems or condensed phases non-additivity can become substantial [8]. A special example of system involving complicated interplay between many forces of different character [START_REF] Saenger | Principles of Nucleic Acid Structures[END_REF][10][11] is found in the case of stabilization of polynucleotide chains. Many contributions are important for DNA involving contacts of its constituents as nucleobases, sugar moieties, phosphate backbones, surrounding water molecules, counter ions and diversity of ligands potentially present in solutions at physiological conditions. Even interactions between two heterocyclic nucleobases are non-homogeneous and typically are divided into three classes, namely hydrogen bonding, intra-and inter-strand stacking. It has been suggested that interaction between nucleic acids bases can in some extend exhibit such nonadditive nature [12][13][14][15][16]. Although base-base interactions in dinucleotide steps of B-DNA are usually pair-additive but in case of d(GpG) dinucleotide steps the absolute value of manybody contribution can be large and significantly structure-related [12][13][14]. As it was documented previously [12][13][14][15][16] high non-additive character of intermolecular interaction energy (IIE) has been found for many d(GpG) dinucleotide steps in B-DNA conformations.

Although three-and four-body terms posses opposite signs they do not cancel each other and in most cases positive non-additivity is observed [15,16]. The origin of such unique properties of d(GpG) steps is expected [12][13][14][15][16] to come from electrostatic repulsion between polar guanine-cytosine pairs what also was related to polymorphic changes of guanine-rich sequences of DNA double helices [16][17][18][19]. On the other hand much weaker electrostatics occurring in d(ApT) steps do not introduces significant non-additivity. However, small values of non-additive terms were also noticed [13,14] for steps comprising two AT pairs. Many molecular clusters comprising for example water, ammonia or formamide molecules [2][3][4] exhibit similar sensitivity to polarization effects. Unfortunately, the level of theory used for characteristics of intermolecular interactions is a crucial factor determining magnitude of binding energies. This is related not only to level of electron correlation estimation [11][12][13][14][15][16][17][18][19], corrections of basis superposition errors (BSSE) [20] and quality of applied basis set [11][12][13][14] (or its extrapolation to complete one) but also to non-additive contributions to molecular cluster stabilizations [12][13][14][15][16]. As it was nowadays documented [16] accurate values for manybody effects require extensions of electron correlation beyond MP2 levels since this approximation can only encounter non-additivity via correlated exchange and induction components [21,22]. However, due to fortunate compensation of errors such levels of theory as MP3, L-CCD or CCSD lead to substantially the same predictions as CCSD(T) one [16].

Due to the computer resources required for characteristics of nucleobases clusters there are only few papers quantifying the non-additive nature of intermolecular interactions [13][14][15][23][24][25]. Unfortunately the size of analyzed systems prohibit application of advanced levels of theory as routine computations for broad range of nucleobases conformations. Bearing in mind formal restrictions of MP2 level [21,22] the qualitative estimation of non-additive nature of dinucleotide steps was characterized here for all ten unique dinucleotide steps adopting B-DNA conformations. Let consider schematic representation of model dinucleotide step as it is presented in Fig. 1. The total binding of four nucleobases in such system can be written in the following general form:

4 ( ) ( ) 1
total tetr tetr d XpY d XpY i i E E E      (1)
The 5'-X and 3'-Y monomers present in the strand I are numbered as 1 or 2, whereas corresponding complementary bases placed in the strand II are labeled as 4 or 3, respectively.

From the formal point of view all terms are to be estimated in the basis sets of the whole tetramer, what was indicated by tetr superscript. The decomposition of the total interaction energy into clusters of multi-body expansion can be defined as follows [26,27]:

44 (2) (3) (4) ( ) ( ) 11 total d XpY ij ijk d XpY j i k j i E              (2)
The analyzed dinucleotide steps are characterized by the following six two-body terms describing pairwise interactions between monomers:

(2) {i,j i}={1,2,3,4}

tetr tetr tetr ij ij i j EEE      (3)
The subscripts denote monomers according to notation used in Fig. 1. For the dinucleotide steps there are four tree-body terms defined as follows:

(3) (2) (2) (2)
{i,j i,k j}={1,2,3,4}

tetr tetr tetr tetr ijk ijk i j k ij ik jk EEEE             (4)
Finally, the only four body term adopts one straightforward form:

4 4 4 (4) (2) (3) ( ) ( ) 1 1 1 tetr tetr d XpY d XpY i ij ijk i j i k j i EE                 (5)
The additive part of the intermolecular interaction energy is usually defined as the sum of all two-body terms but all higher order components constitute so-called non-additive part. The above scheme, although very simple and straightforward, is rather time-consuming and usually not-applicable for large number of dinucleotide conformations and advanced quantum chemistry levels. Thus, the most obvious simplification is neglecting of all non-additive terms and estimation of all energies exclusively in the basis sets of dimers:

4 ,dim () 1
()

add dim dim dim d XpY ij i j ji E E E E       (6)
In such procedure monomers energies are recomputed separately for each pair but the overall CPU cost is about one order lower than full computations in tetramer basis set. The following difference quantifies the error of pairwise additivity of binding energy:

,dim ( ) ( ) ( ) non add total add d XpY d XpY d XpY E E E       (7)
The aim of this study is the quantitative analysis of magnitude of non-additive contributions to the total binding energy of dinucleotide steps and validation of assumption of pair-wise additivities of dimers stabilization energies. Besides, non-standard dinucleotide steps were considered that comprise one 8-oxoguanine molecule denoted by Z in two alternative sequences d(ZpG) or d(GpZ).

Methods

Model d(XpY) steps were prepared according to procedure applied in our previous investigations [28][29][30][31][32] used for description of energetic heterogeneities in canonical and oxidized central guanine triad of B-DNA telomeric fragments [28,29], polymorphism-related heterogeneities of guanine stacking in B-and A-DNA forms [30], characteristics of inter-and intra-strand stacking interactions in d(CpG) and d(GpC) steps found in B-DNA, A-DNA and Z-DNA crystals [31], quantification of all possible intra-strand stacking interactions between nucleobases [32] and assessment of many-body contributions to d(GpG) and d(CpC)

dinucleotide steps [15]. The essential feature of applied approach is the characteristics of statistically significant number of conformations leading to distributions of intermolecular interaction energies (IIE) along with standard deviations and other statistical parameters.

According to Standard Reference Frame [33] the complete description of three-dimensional arrangements of bases in nucleic acid structures is achieved by definition of values of base pair (shear, stretch, stagger, buckle, propeller and opening) and base pair step (shift, slide, rise, tilt, roll and twist) parameters. The values of these parameters corresponding to B-DNA crystals were taken from the Nucleic Acid Database [34] for all dinucleotide steps comprising canonical nucleobases. Preparation of model dinucleotide steps was done with an aid of X3DNA program [35], in which fragments library was modified by replacement of native geometries (without hydrogen atoms) with monomers pre-optimized on MP2/aug-cc-pvDZ (aDZ) level (imposing C s symmetry). These monomers along with sets of 18 parameters (six for either hydrogen bonded pair and additional six for stacked two-bases pairs) were used for preparation of tetramers in which orientations of nucleobases exactly correspond to ones found in the original B-DNA structures. In all cases the whole sugar-phosphate backbones were simplified just by hydrogen atoms. All ten unique model d(XpY) dinucleotide steps were used for single point energy estimation using density fitting MP2/aug-cc-pVDZ (aDZ) level of theory. The counterpoise correction for BSSE error [20] was included in all single point calculations of dimers, trimers and tetramers. The MolPro package [START_REF] Werner | a package of ab initio programs designed[END_REF] was used in quantum chemistry calculations. Many body contributions were characterized according to formulas presented in the proceeding section.

Results and discussion

According to procedure described above there were analyzed 162 model dinucleotide steps S1 in supplementary material. As it was schematically presented in Fig. 1 there are three distinct types of intermolecular interactions that stabilize model dinucleotide steps. The most dominant contribution to the total binding comes from two-body interactions typically described as hydrogen bonding, inter-and intra-strand stacking. Three-body interactions consist of four distinct terms involving all combinations of intra-strand stacked pair interacting with one nucleobase via hydrogen bonding. Finally, there is only one fourbody term.

The magnitude of many body contributions

The calculated values of non-additive contributions to intermolecular interaction energy in ten possible model dinucleotide steps are presented in Fig. 2 and3. First of all it is visible distinct nature of d(GpG) steps since in most cases significant non-additive part is obtained. For this step the distributions of both three-and four-body terms can vary significantly as it is presented in Fig. 2. The strongest attraction coming from three-body terms has been found for d(GpG) dinucleotide step in B-DNA coded as BD0085 and it is as high as -18.9 kcal mol -1 .

The corresponding value of four-body term is also unusually high and equals +20.2 kcal mol - 1 . On the contrary in another utmost case (BD0015) it is noticed extremely high repulsion coming from three-body term (+24.2 kcal mol -1 ) and associated attraction originating from four-body term (-21.2 kcal mol -1 ). However, these two examples are rare extremes since more than 90% of analyzed d(GpG) steps are characterized by three-body terms lower than +6 kcal mol -1 . The corresponding span of four-body contributions is from -5 kcal mol -1 to +2 kcal mol -1 . Closer inspection [15] reveals that not all three-body contributions are of the same magnitude for d(GpG) steps. All terms, except 134 one are usually quite small and typically do not exceed 2.0 kcal mol -1 . However, interactions stabilizing 134 trimer have completely different nature. These contributions can have magnitude of several kcal mol -1 and are the main source of non-additivity in d(GpG) steps. This three-body cluster involves stacking of two cytosine molecules in strand II and hydrogen bonding of 5'-guanine present in strand I.

Interestingly, although three-body contributions may have attractive or repulsive character, always opposite sign of four-body term is observed. However, this does not lead to cancelation of many body contributions and in case of d(GpG) step the total non-additivity is repulsive. As it was presented in Fig. 2

Linear relationships of intermolecular interactions

Since the intermolecular interactions between canonical and modified nucleobases are important aspect characterizing polynucleotide chains there is great temptation in assuming pairwise additivity of interactions stabilizing tetramers. This is justified not only by economy of computations but also by results presented in the previous section. However, in some cases, as for example for extreme conformations of d(GpG) steps, non-additive contributions can be significant. Thus, it is interesting to analyze the relationship between exact values of dinucleotide stabilization energies and simplified ones obtained after neglecting of all non-additive terms in the binding energy of d(XpY) steps and utilization of exclusively dimer basis sets. Thus, the following formula is proposed for regression analysis:

, ,dim ( ) ( ) total est add d XpY d XpY E a E b     (8)
where a and b are regression coefficients and left-hand side estimates the total binding of nucleobases in d(XpY) step as the best fitted values to exact IIE of sum of pairwise additive terms estimated in dimer basis set. The applied regression analysis leads to plots provided in Fig. 4. As one could expect there is observed quite acceptable linear correlation between total binding energy and one estimated as simple sum of pairs energies (in the basis set of dimers).

Even extreme cases for which dinucleotide steps are energetically unfavorable (+3.28 kcal mol -1 in case of d(ApC) found in BDL035 or with very small attraction -6.1 kcal mol -1 for d(ApA) (present in BDJ081) fulfill linear-dependence. At first glance one could recommend utilization of only one linear regression function in Fig. 4, but previously analyzed mean values of all non additive terms (see Fig. 3) suggested that three different classes of dinucleotide steps are to be considered. Since regressions presented in Fig. 4 allow for estimation of not only the total binding but also non-additive part it is interesting to see if linear relationships also hold for the latter quantity. Interestingly, due to different origin of non-additive parts in dinucleotide steps application of only one regression line would lead to very poor estimation of non-additive terms. On the contrary, as it was presented in Fig. 5, three relationships lead to quite acceptable description of non-additive terms by exclusively pairs interactions. Hence, there are three classes of dinucleotide steps with respect of manybody contributions to the total binding energy in dinucleotide steps. This can be further justified by values presented in

Additivity of interactions in dinucleotide steps comprising 8-oxoguanine

Many products of nucleobases oxidation were identified in numerous irradiation and chemical experiments both in vivo and in vitro conditions of DNA and DNA-protein complexes (chromatin) [START_REF] Von Sonntag | Free-Radical-Induced DNA Damage and Its Repair, A Chemical Perspective[END_REF][START_REF] Rice-Evans | Free Radicals and Oxidative Stress : Environment, Drugs and Food Additives[END_REF]. The use of gas chromatography-mass spectrometry with selected-ionmonitoring technique proved steady state levels of many free-radical derived analogues of canonical nucleobases [START_REF] Halliwell | Free Radicals in Biology and Medicine[END_REF]. They are believed to play crucial role in ageing and different stages of cancer development [START_REF] Von Sonntag | Free-Radical-Induced DNA Damage and Its Repair, A Chemical Perspective[END_REF][START_REF] Rice-Evans | Free Radicals and Oxidative Stress : Environment, Drugs and Food Additives[END_REF]. The structural and energetic consequences of B-DNA oxidation were not described till now. It is reasonable to expect that presence of 8-oxoguanine in B-DNA can significantly alter intermolecular interactions [START_REF] Cysewski | [END_REF]. For example intra-strand stacking interactions of 8-oxoG with guanine as well as hydrogen bonding with cytosine are stronger compared to canonical guanine at least in telomere repeat units [28,29]. It is interesting to see if modification of guanine leads also to any non-additivities in base-base interactions occurring in oxidized dinucleotide steps. Since there are not available sufficient amount of experimental data comprising 8-oxoguanine alternatively molecular dynamics simulations were used for collecting structural parameters defining oxidized dinucleotide steps. Details are provided elsewhere [START_REF] Cysewski | [END_REF]. The decomposition of dinucleotide stabilizing energies into many-body contributions can be found in Fig. 6. Presence of 8-oxoguanine in d(ZpG) or d(GpZ) dinucleotide steps introduces significant changes upon non-additivity. First of all, unlike to canonical d(GpG) sequence, both three-and four-body contributions are usually of the same sign. Thus, although actual values of non-additive terms are smaller if compared to canonical guanine dinucleotides it is expected that cumulative effect occurs. The mean values of all three-body terms is about +2.0 and +2.5 kcal mol -1 for d(ZpG) and d(GpZ) dinucleotide steps, respectively. The mean values of four body contributions are much smaller and do not exceed 0.13 kcal mol -1 in either case. Thus, many-body contributions are of repulsive character reducing total stabilization of the two stacked pairs. The mean values of non-additivities characterizing oxidized dinucleotide steps can reach 5% of total IIE. This is similar amount to one found in case of d(GpG) steps [15]. No sequence dependence of many body contributions is observed for oxidized dinucleotides. There is another interesting features of three body terms found for modified dinucleotide steps. Majority of 123 and 134 contributions are systematically smaller than 124 and 234 ones. The common patterns in these terms are formed by purine-pyrimidine inter-strand stacking (G 1 |C 3 , Z 1 |C 3 , C 2 |G 4 or C 2 |Z 4 ). It seems that 13 interactions introduces attractions in analyzed steps, while 24 inter-strand interactions are responsible for repulsive contributions irrespectively of the analyzed sequences. In case of d(GpG) this feature was not observed [15] since major contributions to the total non-additivity came from 134 term. Finally, it is worth mentioning that there are linear relationships between full estimation of IIE values and ones utilizing pair-wise additive assumption. The applied linear regression analysis led to quite acceptable correlation with R 2 =0.983 for oxidized dinucleotide steps and R 2 =0.996 for canonical ones. This relationship is presented in subsequent study of this issue [START_REF] Cysewski | [END_REF]. Finally, the total binding of four nucleobases in oxidized d(GpG) steps is about 2 kcal mol -1 stronger than found for canonical dinucleotide.

Conclusions

Results presented above qualitatively characterize non-additivity of base-base interactions stabilizing d(XpY) dinucleotide steps. Although applied model excluded sugar-phosphate backbone, the most important features of intra-and inter-strand stacking along with hydrogen bonding were properly preserved. The most important aspects of presented data is rational representation of four bases conformations and adequate level of computations. The first facet relies on the representation of meaningful conformations matching to ones occurring in B-DNA crystals. The main conclusion drawn from presented results is that pair-wise additivity simplification is rational since corresponding intermolecular interaction energy of d(XpY) steps suffers far less than 5% if compared to correct values. The four-and three-body terms are small and usually of the opposite signs, what further reduces its impact on the total binding energy. Dinucleotide steps comprising two stacked guanine molecules seem to be the only canonical systems for which non-additivity is non-negligible. However, as it was demonstrated in the previous report [15], in this case there are observed linear relationships which allow for corrections of IIE values without necessity of energy computations for the whole tetramer. Despite the fact that many body contributions are small, they can be found as interesting source of classification of dinucleotide steps into three distinct classes. Hence, despite analogical distributions of energetic asymmetry the high non-additive contributions are observed only for canonical and oxidized d(GpG) steps. This suggest that the selection of dinucleotide conformations is adequate and presented conclusions are expected to be consistent also for other conformations of dinucleotide steps. The results presented in this paper are preliminary and further investigations are necessary. Several aspects were excluded from the analysis that probably can affect both stabilization energies of dinucleotide steps and many-body contributions. Among them the applied level of computations is the first way of improvements. More sophisticated methods should be applied for proper inclusion of electron correlations. Although, presented results embrace only small portion of this effect due to applied level of computations it is reasonable to expect that extension of the method will not change the general conclusion of the magnitude of manybody contributions do dinucleotide steps energetic. Besides there are other aspect worth additional considerations. For example it is commonly known that nucleobases comprising amino groups are not planar. These exocyclic substituent exhibit pyramidalization, which probably affects both IIE and many-body terms. Also relaxation of monomers geometry, especially hydrogen bonds are worth of precise evaluation. The importance of all these aspects is to be validated in further investigations.

Tables

Figure captions

add dim dim dim d XpY ij i j ji E E E E        ) Fig. 5
The applicability of linear relationships obtained in Fig. 4 for estimation of nonadditive contributions based on pair interactions estimated using dimer basis sets The distribution of energy asymmetry coefficient (δ) estimated for all of analyzed dinucleotide steps

7 )

 7 Page 4, end of the first paragraph "The mean values of all ..." should be "The median values of all ...". Thank you! 8) Page 4, end of the first paragraph I don't think the following statement is true: "On the contrary for such steps as d(ApA), d(TpA) and d(ApT) the nonadditivty has attractive nature." From Figure 2 it follows, that only AA is attractive, while AT and TA are repulsive. Similarly, the third class does not have all steps with repulsive nonadditivy (AA has attractive nonadditivity) I have also updated Fig.3 and changed the text accordingly. 9) Concerning the third class, small value of nonadditivy term in steps with two AT pairs was already noticed e.g. in Ref 11, or recently in DOI: 10.1021/jp910788e. I think this should be acknowledged by proper reference.

  including 48 of d(GpG), 23 of d(ApA), 16 of d(CpG), 15 of d(TpA), 10 of either d(GpC), d(GpA), d(ApG), d(ApC), d(CpA) and d(ApT). All analyzed dinucleotide steps were described in Table

  many-body contributions found for all other dinucleotide steps are much smaller and can be assumed as negligible since three-body terms are always within -1.8÷1.7 kcal mol -1 and four-body ones belong to interval -0.7÷0.2 kcal mol -1 . Since stabilization energy of model dinucleotide steps usually exceed -40 kcal mol -1 such small non-additivity affects at most modestly the overall biding of four nucleobases in d(XpY) steps. The median values of all non-additive contributions were presented in Fig. 3. Although all medians are positive the dinucleotide steps can be divided into three classes. As it was mentioned above the d(GpG) steps are characterized by highest repulsion coming from non-additive interactions. On the contrary for d(ApA), d(TpA) and d(ApT) dinucleotides nonadditivity although positive is close to zero. The third class comprises the rest of dinucleotide steps, for which non-additivity is quite modest and repulsive. The statistical analysis applied to these three classes of dinucleotide steps proved that the differences in the median values are greater than would be expected by chance (P = <0.001). Hence all three dinucleotide distributions are statistically different. The normality test failed only for the first class comprising d(GpG) dinucleotides. The other two sets passed both normality (P = 0.134) and equal variance tests (P = 0.096).

  polarity. Thus, polarity expressed in terms of dipole moments along with electron correlation and SCF contribution to the total binding energy are quite satisfactory criterions of classification of dinucleotide steps into three distinct classes. Although the absolute values of non-additive terms seem to be negligible they are still interesting source of distinguishing of dinucleotide steps.

  First class comprised d(GpG) steps that were found as most polar among all dinucleotide steps. To the second class one can include d(ApA), d(TpA and d(ApT) steps. They are characterized by highest values of SCF contributions to the total IIE and lowest polarities. The rest of dinucleotide steps constitute the third class with modest polarity and electrostatic contribution to the total interaction energy. Oxidation of guanine in d(GpG) dinucleotide steps introduces non-negligible non-additivity to IIE that is of the same order as observed for canonical systems. In the case of many potential structures contribution to analyzed quantity there is always question related to arbitral selection of analyzed systems. It is impractical to analyze hundreds of d(XpY) conformations and misfortune sampling may significantly affect conclusions. However, in the previous paper[15] it was demonstrated that there is simple origin of non-additivity in the case of d(GpG) dinucleotide steps. The interactions between hydrogen bonded nucleobases as well as inter-strand stacking interactions can vary significantly for different conformations of four bases. As the quantitative measure of such feature the energy asymmetry coefficient (δ) is proposed. It is defined as the sum of two following terms. First one denotes the difference in hydrogen bonding energies between two neighboring complementary nucleobases and the second one characterizes difference between 5-G|C'-5' and 3-G|C'-3' inter-strand stacking interactions. In case of d(GpG) dinucletode steps this parameters emphasizes the importance of hetero-bases contacts. There is straightforward way of extension of δ coefficient for other tetramers and characterize the symmetry/asymmetry diversities (in energy sense) of all dinucleotide steps. In any case close to zero values of δ coefficient indicate high symmetry of intermolecular interactions. On the contrary the positive values of this parameter suggest stronger interactions on 3' side and vice versa negative values emphasizes stronger attraction on 5'-side. In Fig.7there is presented distribution of δ values estimated for analyzed dinucleotide steps. As is it visible all analyzed systems may adopt wide range of energetic asymmetry covering broad range of structures.

Fig. 1 Fig. 2 Fig. 3 Fig. 4

 1234 Fig. 1The schematic representation of intermolecular interactions in model d(XpY) dinucleotide steps. There are three distinct types of contacts, namely intra-strand stacking (5'-X 1 /Y 2 -3', 5'-Y' 3 /X' 4 -3'), inter-strand stacking (5'-X 1 |Y' 3 -5', 3'-Y 2 |X' 4 -3') and hydrogen bonding (5'-X 1 -X' 4 -3', 3'-Y 2 -Y' 3 -5'), where X,Y denotes one of four nucleobases and X',Y' stands for corresponding complementary base

Fig. 6

 6 Fig. 6Distribution of many body contributions to intermolecular interaction in oxidized d(GpG) dinucleotide steps. Open circles denote d(GpZ) steps, while black circles stand for d(ZpG) sequence, where Z denotes 8-oxoguanine

Figure 6 .

 6 Figure 1 Click here to download high resolution image

Figure 6

 6 Figure 6 Click here to download line figure: Fig6.doc

2 Figure 1 .

 21 Figure 7 Click here to download high resolution image

Figures 1 Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 12345 Figures 1-7 (alternative file format) Click here to download attachment to manuscript: JMM01145R1_Figures.doc

Table 1

 1

	, which collects the mean values of SCF and electron
	correlation contributions to the total values of IIE. Besides, the mean values of dipole
	moments of analyzed dinucleotide steps are also provided. Interestingly, although all SCF
	energies are negative indicating strong electrostatic contribution to the total binding energy,
	there are observed quite significant diversities between particular dinucleotide steps. First of
	all d(ApA), d(ApT) and d(TpA) systems are characterized by highest values of SCF term.
	Additionally, these tetramers belong to the less polar steps. On the contrary the d(GpG) is
	known [11,12] as the most polar system among all ten possible dinucleotide steps, what also

is confirmed by mean values of dipole moment collected in

Table 1. Besides, much lower values of SCF contribution are noticed in this case if compared to previously described set of dinucleotides. The rest tetramers constituting class of dinucleotide steps that is characterized by much stronger attraction coming from electron correlation energy and are of modest

Table 1

 1 The mean values of dipole moments, SCF and electron correlation contributions to the total binding energy of model dinucleotide steps. Energies are expressed in kcal mol -1 and dipole moment are given in Debyes. Bold face was used for demonstration of differences

	between three distinct classes of dinucleotide steps
	class	step	E 	, () total SCF d XpY	E 	, () total correl d XpY	dipole moment
	I	d(GpG)	-16.9±13.1	-51.0±10.7	11.8±1.3
	II	d(CpG)	-24.7±16.0	-56.7±13.6	2.7±1.3
		d(GpC)	-29.2±1.9	-64.9±1.3	4.8±0.4
		d(GpA)	-13.2±6.9	-47.7±4.7	7.2±0.3
		d(ApG)	-11.3±13.5	-44.2±12.0	8.4±0.3
		d(ApC)	-6.2±23.5	-44.2±16.8	4.5±0.5
		d(CpA)	-23.3±2.2	-52.0±0.4	4.5±0.3
	III	d(ApA)	-2.2±6.3	-38.3±4.7	3.9±0.4
		d(ApT)	-0.9±4.7	-37.6±5.6	1.6±1.3
		d(TpA)	0.3±12.2	-32.3±8.4	1.7±0.4

Acknowledgments

Results were obtained the part of computational grant no 39 of Poznań Supercomputing and Networking Center (PSNC, Poznań, Poland). The allocation of computational facilities are greatly appreciated.

Supplementary materials

Table S1

The list of pdb files deposited in Nucleic Acid Database [33,34] that were used in this study as the source of structural parameters of dinucleotide steps. First two columns comprise B-DNA codes according to NDB and PDB identification styles. Resolution is provided in Angstroms. The final column comprises dinucleotide acronym along with the position of suitable base pairs in the B-DNA chain. If alternative pdb structures were used with to the same ID it is indicated in the bracket.

NDB


J Riegler
email: j.riegler@ucl.ac.uk

B Allain

R J Cook

M F Lythgoe

Q A Pankhurst

Magnetically assisted delivery of cells using a magnetic resonance imaging system

A simple analytical model is presented which enables rapid interactive prediction and control of magnetically labelled cells in an arterial bifurcation using magnetic field gradients produced by a magnetic resonance imaging system (MRI). This model is compared against experimental results for human mononuclear cells labelled with micron sized superparamagnetic iron oxide particles. Experimental and theoretical results highlight the importance of cell aggregation for magnetic targeting in a strong magnetic field. These predicted aggregates are confirmed via confocal endoscopy which allows the visualisation of cell aggregates and their movement inside a vascular flow model in a 9.4T preclinical MRI scanner.

1.

Introduction

The increased understanding of cellular processes and the role of cells during embryonic development and regeneration has opened new research avenues for cell transplantation therapies. These cell therapies use a variety of cell types from terminally differentiated cells such as T-cells [START_REF] Nakajima | A phase I study of adoptive immunotherapy for recurrent non-small-cell lung cancer patients with autologous gammadelta T cells[END_REF] to embryonic stem cells [START_REF] Corporation | World's first clinical trial of human embryonic stem cell therapy cleared Regen[END_REF]. Several clinical trials are already investigating the potential of adult stem cells contained in the bone marrow for cardiac repair after myocardial infarction as well as other possible applications [START_REF] Brinkman | Autologous stem cell transplantation in children with severe progressive systemic or polyarticular juvenile idiopathic arthritis: Longterm followup of a prospective clinical trial[END_REF][START_REF] Snowden | Clinical stem cell therapies for severe autoimmune diseases[END_REF][START_REF] Yeo | Autologous bone marrow-derived stem cells for ischemic heart failure: REGENERATE-IHD trial Regen[END_REF][START_REF] Yerebakan | Autologous bone marrow stem cell therapy for the ischemic myocardium during coronary artery bypass grafting Minim[END_REF][START_REF] Wollert | Intracoronary autologous bone-marrow cell transfer after myocardial infarction: the BOOST randomised controlled clinical trial[END_REF].

The efficacy of such treatments is determined by many factors such as cell type, differentiation potential, delivery route, delivery time, local micro-environment, inflammation, cell retention etc. Current approaches for cell delivery and retention frequently end up achieving less than 10% retention of administered cells [START_REF] Tang | Intracoronary administration of cardiac progenitor cells alleviates left ventricular dysfunction in rats with a 30-day-old infarction[END_REF]. Magnetically assisted delivery is one of a range of strategies proposed to address this issue.

Magnetic targeting strategies rely on the fact that cells which have been made responsive to magnetic fields by internalisation of iron oxide particles can be steered over a distance using magnetic fields [START_REF] Pankhurst | Applications of magnetic nanoparticles in biomedicine[END_REF]. Particles used for such applications exhibit superparamagnetism due to their size. As other cells in the body do not contain sufficient amounts of magnetic materials only labelled cells are affected and can be concentrated and retained in the target area. In order to concentrate magnetic particles or cells in an area, a magnetic field and a field gradient is necessary [START_REF] Pankhurst | Applications of magnetic nanoparticles in biomedicine[END_REF]. These can be generated by a number of sources including permanent magnets outside of the body, implanted magnets or permanent/electro magnets external to the body in combination with ferromagnetic materials such as needles or stents to generate high field gradients. A general limitation of magnetic targeting, especially for permanent magnets and electromagnets positioned outside of the body is the rapid decay of the magnetic force with increasing distance from its source, which follows an inverse square relation. This rapid decay of the magnetic field strength and gradient strength, leads to a parallel decay of magnetic forces, limiting its applicability to organs close to the body surface [START_REF] Chen | Capture of magnetic carriers within large arteries using external magnetic fields[END_REF][START_REF] Kim | Targeting to carcinoma cells with chitosanand starch-coated magnetic nanoparticles for magnetic hyperthermia[END_REF][START_REF] Hofmann | Combined targeting of lentiviral vectors and positioning of transduced cells by magnetic nanoparticles[END_REF][START_REF] Kyrtatos | Magnetic tagging increases delivery of circulating progenitors in vascular injury[END_REF][START_REF] Lübbe | Physiological aspects in magnetic drug-targeting[END_REF].

A potential alternative was proposed and explored by Mathieu and Martel [START_REF] Martel | Automatic navigation of an untethered device in the artery of a living animal using a conventional clinical magnetic resonance imaging system[END_REF][START_REF] Martel | Interventional procedure based on nanorobots propelled and steered by flagellated magnetotactic bacteria for direct targeting of tumors in the human body[END_REF][START_REF] Mathieu | MRI System as a mean of propulsion for a microdevice in blood vessels[END_REF]. In experiments on live pigs, they were able to actuate and steer an intravenously introduced 1.5 mm steel ball through parts of the porcine vascular network by manipulating the gradient coil currents of a standard 1.5T magnetic resonance imaging (MRI) system. By doing so they were able to successfully demonstrate the concept of magnetic resonance targeting (MRT), viz. magnetically assisted delivery using an MRI system; although it could be argued that they did not establish the feasibility of MRT in man, as the use of such large steel objects is inconceivable due to the risk of thrombosis.

Nevertheless, if a variant of MRT could be made to work in man, it could be highly advantageous and useful. To this end, we have been exploring the use of magnetically labelled cells. In recent work we have shown that substantial actuation forces can be generated on such cells [START_REF] Riegler | Targeted magnetic delivery and tracking of cells using a magnetic resonance imaging system[END_REF], at a level that makes them suitable for MRT. In principle one can then imagine a viable human therapeutic procedure in which, for example, a patient's own stem cells could be harvested, magnetically labelled, and then re-introduced into the bloodstream, where their unchanged membrane characteristics would ensure biocompatibility, while the internalised magnetic particles would render them suitable for magnetically assisted delivery. Furthermore, there are at least four major advantages in using magnetically labelled cells and the gradient coils in MRI systems for targeting applications. (1) The high magnetic field strength (more than 1.0 T) ensures the full magnetisation of the magnetic particles in the cells, and maximises their actuation potential. [START_REF] Corporation | World's first clinical trial of human embryonic stem cell therapy cleared Regen[END_REF] The high spatial and temporal control over the gradient strength and direction in MRI systems allows high-precision and active guidance of particles along a chosen path. (3) The conventional use of MRI systems to generate high contrast soft tissue images can be combined with the targeting application to both guide the magnetically labelled cells to the area of interest and to monitor their accumulation at the target over time. (4) MRI systems are widely available in hospitals and clinics all over the world, so that the development of a delivery method based on MRT would be similarly widely available, without the need for major capital expenditure on new equipment.

We are therefore continuing our exploration of the feasibility of MRT on magnetically labelled cells as a potential human therapeutic intervention. In the current study we focus on an important aspect of the magnetically assisted delivery process itself, namely the ability to understand and predict the direction of motion of magnetically labelled cells in the vasculature. We start from the principle that the best way to achieve this is to derive a simple analytical model, and preferably an algebraic one, to allow rapid and interactive processing as part of the MRT process. In the following we show that this is possible, and we derive a simple analytical model describing magnetic cell steering in an MRI scanner. We introduce several practical concepts to assist in the cell tracking and steering process, including a weighting factor to account for both the cell distribution in a tube that is projected into a 2D plane, and the presence of prolate spheroids to account for cell aggregation. We compare this model against experimental results and highlight the importance of cell aggregation at lower flow rates. In order to confirm these predictions we used confocal endoscopic imaging inside and outside the bore of the MRI scanner to visualise cell distribution and aggregation.

2.

An analytical model of magnetically assisted cell delivery in an MRI system

The movement of suspended cells in a magnetic field is primarily governed by the magnetic force, drag force, gravity, inertia and magnetic dipole interactions. The model described below takes the magnetic force and drag force into account while inertia and gravity are ignored. Magnetic dipole interactions are indirectly included by assuming cells aggregating to bigger objects. The assumption of cell aggregation and factors influencing it, such as cell concentration, viscosity and fluid velocity will be experimentally parameterised and assessed later (see section 3.1 and 3.2), while the influence of aggregate shapes will be theoretically assessed (see section 4.5).

2.1

Magnetic field The main driving force for cell movement inside of the MRI scanner is the field gradient produced by the coils used to generate the imaging gradients. The magnetic field can be described using Maxwell's equations:

, and

(2) where µ and are the magnetic and electric permeabilities of free space, the electric field and the current density. The bore of an MRI scanner is a source free field which means that there are no electric fields or currents present. Accordingly equation 2 reduced to: .

(

) 3
Equations 1 and 3 can be rewritten in their component forms as:

, and

Equations 4 and 5 contain five independent derivatives. The first three of them can be selected as the linear field gradients which are derivatives of as the main field is orientated along the z-axis in MR systems:

! "# $ " ; % ! "# $ "% ; ! "# $ " . (6)
For a symmetrical z gradient ( ) in a cylindrical coordinate system ' ( ') and ' % ( '* should be identical. We can therefore introduce a symmetry parameter :

! "# + (" , $ ; - ! "# . ("% , $ . ( 7
)
As the last independent variable we can chose equation 5c:

! "# . " "# + "% . ( 8
)
All the partial derivatives described in Equations 4 and 5 can be expressed with the independent variables described in Equations 6, 7 and 8. The field produced by the gradient coils can then be expressed as:

/ % 0 1 2 2 2 3 "# + " "# + "% "# + " "# . " "# . "% "# . " "# $ " "# $ "% "# $
For MRI gradient coils < and due to the cylindrical geometry -(=. Equation 9 can therefore be simplified to:

/ % 0 > ? @ ) 8 ? @ * % 8 ) % * 8 A , (10)
which is an especially convenient algebraic expression for the magnetic field in MRI systems [START_REF] Han | Gradient waveform synthesis for magnetic propulsion using MRI gradient coils[END_REF], as it is common practice for scanner manufacturers to supply data on , % , and amongst the specifications of their instruments (for field maps see figure 1).

In addition Equation 10shows explicitly that gradients along the x-and y-axes are produced by the field. Hence we can have a magnetic force along those axes despite the magnetisation of material introduced into the field of the scanner being solely along the z-axis.

2

Magnetic force on cells The magnetic force acting on a particle is dependent on its magnetic susceptibility, the field strength and the gradients of the field, which can be expressed as:

B C :D ; , and (11)
D EFG H I ( 12
)
where J is the particle volume and FK is the magnetic susceptibility difference between the particle and its surrounding fluid. This equation can be simplified for the case of MRT of cells since the magnetic particles will be fully magnetised due to the strong L field of the scanner. We can therefore simplify equations 11 and 12 to:

B C MJ:N O P; (13)
where M is the number of magnetic particles in a cell and N O is the saturation magnetisation for a magnetic particle. As the magnetisation of particles in an MRI scanner is mainly along the z-axis, the x and y magnetisation components can be neglected reducing Eq. 13 to:

QR CS R CS R CS T U MJV W X U . (14)
For our experiments only gradients across the x-axis were used ( ). We can therefore reduce our mathematical model to two dimensions and simplify Eq. 14. An additional scalar needs to be added to account for the fact that imaging gradients can only be switched on for short times at high amplitudes and are therefore operated in an on/off fashion, referred to as the systems duty cycle (Y Z ). Using that and Equations 5b we can write the magnetic force along the x-direction as:

R CS MJY Z V [ . ( 15
)
The magnetic force is opposed by the drag force of the fluid which can be stated as Stokes drag force for low Reynolds numbers:

R \ ]^_`Zab Z b c d , (16)
where `Z is the cell radius, _ the dynamic viscosity and b Z b c is the difference between the cell velocity and the fluid velocity.

Cells are likely to aggregate in chain like structures. We can introduce a form factor e f to extend Stokes drag to non spherical objects [START_REF] Leith | Drag on Nonspherical Objects[END_REF]:

R \ g^_h f ab Z b c de f , and (17)
e f W ? i @ i \ j \ k X . ( 18
)
where h f is the diameter of a sphere with the same projected area as the area of the object projected normal to its direction of motion while h [ is the diameter of a sphere with the same total surface area as the object. For a prolate spheroid Equation 18 can be written as:

e f l ? i @ i m n o pq r stu vw x @mZ y and z {- m n Z n , ( 19
)
where | is the equatorial radius and } the polar radius of a prolate spheroid (a<c) as shown in Figure 1.

We can use Eq. 17 to investigate the effect of cells arranged to prolate spheroids with different diameter to length ratios ~• on targeting efficiency:

~• ! Z m . (20)
As the flow rates used for our experiments are all in the laminar flow region, the inflow profile can be described as:

b c :€; =b c ••• l- ' ƒ @ y , (21)
with ` as the tube radius and b c ••• the average fluid velocity in the tube. The trajectory of cells can be calculated by equating the magnetic and the drag force (see Figure 1). Note that as our model does not include inertial forces, the cells will immediately travel at the maximum velocity determined by the magnetic and drag forces

R C R \ . ( 22
)
The cell velocity in the z-direction is determined by the fluid velocity, as there is no magnetic force acting in the z-direction. Therefore the solution for the particle trajectory along the x-axis will be independent of the solution for z. The substitution of Equations 15 and 17 into Equation 22for x leads to:

\ \ "E… q † j ‡ˆ ‡ i‰Š\ k ‹ k . ( 23
)
The velocity along the z-axis is given by Eq. 21. However as the cell changes its position along the xaxis the velocity changes we therefore have to take the position along x over time into account which leads to:

\ \ =b c ••• l- : ; ' ƒ @ y . ( 24
)
Equations 23 and 24 lead to:

€:OE; € b Z OE . ( 25
)
The velocity of cells due to the magnetic force is given as b Z :

b Z "E… q † j ‡ˆ ‡ i‰Š\ k ‹ k , (26)
and the position along z as a function of time as:

•:OE; =OEb c ••• @ Ž • q n • • •••• i' ƒ n @ n • q • • •••• I ' ƒ n @ • • •••• I n ' ƒ n . ( 27
)
In order to calculate the targeting efficiency, the initial distribution of cells in the tube at the start of the separation length (z=0) needs to be taken into account. It is assumed that the cells are randomly distributed over the circular tube cross-section. If this distribution is projected into the plane we derive:

':€; =` '"" •-' -? ' ƒ , ( 28
)
where ' is a weighting factor for the cell number distribution at the start of the separation area (see Supplementary Figure 1). The trajectories for cells can be calculated using Eqs. [START_REF] Takeda | Development of magnetic force-assisted gene transfer system using biopolymer-coated ferromagnetic nanoparticles[END_REF] and 27 (see Supplementary Figure 2). If they are weighted with Eq. 28 the distribution at the start and end of the separation distance can also be calculated. It is assumed that cells which are above or below the z-axis at the end of the separation length • (see Figure 1) will end up in their corresponding exit tubes. The targeting efficiency can hence be defined at that point as follows:

˜x = ™ " ƒpš›rƒ&ƒoe•r " ƒžƒpŸ ¡ , ( 29
)
where M m•¢x & £¤x is the number of cells at the end of the tube to which the cells were targeted to and M ¥ m¦ is the total amount of cells in both tubes. The subtraction of 0.5 is accounting for the distribution without any steering gradients present, which is 0.5 for both exit tubes.

3

Estimating parameters to achieve a given targeting efficiency In order to estimate how much iron oxide needs to be internalised per cells to achieve a certain steering efficiency Eq. 28 needs to be integrated. For a selected targeting efficiency ˜x, the fraction of M x " ƒpš›rƒ&ƒoe•r " ƒžƒpŸ cells have to end up in the targeted tube. We therefore integrate Eq. 28 from -` to € x if the magnetic force points along the positive x-axis. In this case € x indicates the position along the x-axis for which the integral between it and ` accounts for ˜x of all cells. We can solve the following equation for € x :

M x ` @ ^ ` §€ x {- ƒr n ' ƒ n ` '"" -? ƒr ' ƒ ¨ ' ƒ n ‰ @ . ( 30
)
If we chose ˜x to be 0.9, than form Eq. 28 it follows that M x has to be 0.95. In order to achieve the 0.9 targeting efficiency, a cell starting at € x and z=0 has to reach € at the end of the separation length • . Equation 27 and 25 can be solved for OE:

OE i' ƒ n I @• • ••••ai' ƒ n -ƒr n d . ( 31
)
From Equations 30 and 25, M © can be calculated as well as the field gradient strength or duty cycle necessary to reach a given steering efficiency. For parameters used in the model see Table 1.

3.

Materials and Methods

Cell-steering experiments were conducted with different flow rates, cell concentrations and fluid viscosities to investigate their effects on targeting efficiency. Additionally, confocal endoscopy was performed to investigate parameters used in our model (section 2) describing cell aggregation.

3.1

MRT of labelled cells A vascular bifurcation model consisting of a cylindrical tube with an internal diameter of 0.8 mm and a 30°bifurcation (see Figure 1) was placed into the isocentre of a 9.4T MRI scanner (VNMRS, Varian Inc. Palo Alto, CA) equipped with a 60 mm bore size gradient set (rise time 5 T m -1 ms -1 , max. gradient strength 1 T/m) and connected to an infusion pump (PHP2000 Harvard Instruments, UK). Labelled cells (see Appendix A) suspended in PBS containing 3% foetal bovine serum (FBS) (or 50% FBS) where flushed through the flow phantom infused at flow rates of 0.06 mL/min to 2.5 mL/min, leading to a mean velocity of 0.2 cm/s to 8 cm/s assuming parabolic laminar flow. During the infusion, gradients (amplitude + or -500 mT/m) were applied in the x-direction, perpendicular to the direction of flow (see Figure1). Due to hardware capabilities the gradients were pulsed (2 ms on, 7 ms off using max. slew rate). Cell suspensions leaving each bifurcation outlet tube were collected and cell concentrations were estimated using a haemocytometer. This was repeated three times for each gradient direction and fluid velocity.

3.2

Confocal endoscopy The Cellvizio (Mauna Kea Technologies, Paris, France), a fibered confocal microscope, was used for the acquisition of dynamic microscopic in vitro images of the fluorescent cells. An MR compatible (7 metre long) 600 micron diameter microprobe, comprising approximately 10,000 coherent optical fibres was scanned fibre by fibre via the MKT laser scanning unit behaving as the confocal pinhole and objective lens of a conventional near real-time confocal microscope for remote acquisitions at 12 frames per second. The MR-compatible ProFlex S-microprobe was used which provided an image resolution of approximately 4 µm and a circular field of view with a diameter of 500 µm. The end of the optical fibre was placed into the middle of the inflow channel of the flow phantom where a hole was drilled to connect it to the fluid (see Figure 1 and3). The laser scanning unit emitted a laser beam at an excitation wavelength of 488 nm and collected the fluorescence signal in a bandwidth of 500 nm to 650 nm. Finally, ImageCell Software allowed for real-time visualisation of the microscopic images and for manual settings of the contrast.

4.

Results and discussion

In the following sections we compare our experimental results with model predictions and assess some of the assumptions used in the model, in particular with regard to cell aggregation. Experimental variation of some of these parameters and their effect on targeting efficiency will also be discussed.

Theoretical prediction for single cells and aggregates

Figure 2A shows a comparison between experimental results and theoretical predictions for two different cell concentrations, together with the theoretical prediction without cell aggregation. This figure demonstrates that the theoretical prediction for single cells differs significantly from experimental results, indicating that aggregation could play an important role, which may explain the discrepancy with the experimental data. When cell aggregation was assumed, cells were modelled as prolate spheroids with a polar to equatorial diameter ratio of ~• ¡. Cell aggregation and aggregate shape were included in our model via Equations 15 and 18. The number of cells N in such an aggregate was derived via fitting the experimental data to the model while the diameter length ratio was estimated from confocal endoscopy data. For N with the range 60 to 80 coefficients of determination (R 2 ) were found to be between 0.93 and 0.97, showing a good agreement between the model and experimental data.

It should be stated that the numerical validity of the results is dependent on the validity of the model assumptions of no gravitational and inertial effects. Neither gravity nor inertia are likely to have a significant effect on our predictions, as sedimentation is slow and the cells pass through the flow phantom within a few seconds. Inertia is neglectable for small cells, especially when the magnetic force is applied for the entire time a cell spends in the flow phantom. The assumption of a random distribution for cells entering the separation area is harder to assess as a long tube is used to connect the phantom with the syringe pump. However, applying gradients along the positive or negative x-axis lead to similar results indicating no initial bias. Finally the assumption of a laminar parabolic flow profile is justified for the fluid viscosities and flow rates used in our experiments.

Effect of viscosity and cell concentration

Subsequently we investigated two different cell concentrations to assess if doubling the cell concentration leads to increased cell aggregation and hence increased steering efficiency. This was observed, however the effect was relatively small. The model used a constant aggregate size which is likely to be incorrect particularly for higher flow rates where increased shear forces overcome magnetic dipole interactions and hence aggregation (see confocal results below). We therefore see a discrepancy between the model and experimental data starting between 5 and 6 cm/s for 2x10 6 cells/ml and 6 and 8 cm/s for 4x10 6 cells/ml. Figure 2B shows experimental results and theoretical predictions for two different fluid viscosities. The increases in fluid viscosity lead to a modest decrease in targeting efficiency. Changing the fluid viscosity does affect the drag force which needs to be overcome and to a lesser degree cell aggregation.

This can be seen in figure 2B for which the aggregate size was kept at 60 cells with a diameter to length ratio of 5 and measured fluid viscosities were used to calculate theoretical targeting efficiencies. These predictions fit the experimental data assuming that the aggregate size does not change significantly. The importance of magnetic dipole interactions for targeting applications have been described previously [START_REF] Cregg | Inclusion of magnetic dipole-dipole and hydrodynamic interactions in implant-assisted magnetic drug targeting[END_REF][START_REF] Erica | Particle size, magnetic field, and blood velocity effects on particle retention in magnetic drug targeting[END_REF][START_REF] Mathieu | Aggregation of magnetic microparticles in the context of targeted therapies actuated by a magnetic resonance imaging system[END_REF]. Also there have been previous publications where a similar bifurcation has been modelled but in the context of permanent magnets and without cell aggregation [START_REF] Mohanty | Three-dimensional CFD modelling of a continuous immunomagnetophoretic cell capture in BioMEMs[END_REF][START_REF] Takeda | Development of magnetic force-assisted gene transfer system using biopolymer-coated ferromagnetic nanoparticles[END_REF]. However, in comparison to these studies we used cells with internalised particles and not magnetic particles alone. This changes aggregation behaviour as particles are kept at greater distances due to other cellular compartments and their membranes, while aggregated cells might stick together more readily depending on the cell line used. Forces which oppose aggregation of particles are electrostatic and van der Waals forces. These forces have a shorter range compared to magnetic dipole interactions and are shielded to a large extend by the cells. We should therefore expect slight differences between the aggregation behaviour of ferrofluids and cells even if the same iron concentration per unit suspension volume is used.

Assessment of aggregation using confocal endoscopy

In order to investigate cell aggregation in the context of MRT we used a confocal endoscope. This system operates as a conventional confocal microscope as the coherent fibre bundle is rapidly scanned point by point, conducting light to and from the sample, acting as both the objective lens and the confocal pinhole. Figure 3 shows single frames from videos taken at different fluid velocities without a magnetic field present and with a homogenous magnetic field of 9.4T (bore of the MRI scanner) (see supplementary movie S1-S6). For fluid velocities of 1 cm/s and 2 cm/s (Figure 3A andB) cells travelled as aggregates in the magnetic field of the MRI scanner; these aggregates showed roughly elongated shapes containing between 20-80 cells. We estimated the diameter to length ratio of these aggregates to be around 5. For the same experiments conducted without a magnetic field (outside of the MRI scanner), no cell aggregation was observed. When the fluid velocity was increased to 4 cm/s (Figure3C) cells did not form large aggregates but rather chain-like aggregates with 2 to 10 cells arranged into a chain. The transition between large aggregates (several 10s of cells) and smaller chain like aggregates was observed at a fluid velocity of approximately 3 cm/s. It has to be accepted that the imaging fibre was introduce into the flow phantom (see Figure 3D) and may have caused some minor flow disturbances.

Mathieu et al. [START_REF] Mathieu | Aggregation of magnetic microparticles in the context of targeted therapies actuated by a magnetic resonance imaging system[END_REF] showed that the average chain length of iron oxide particles without fluid flow at a concentration of 0.6 mg/ml and 0.4T was about 31µm and for 1.13 mg/ml at a flow rate of < 1.2 cm/s in an 1.5T MRI scanner approximately 650 µm. In our study, we reached an iron concentration, per unit suspension volume, of approximately 0.5 mg/ml in labelled cells and observe aggregates with a length of 150-250 µm (see S2). This difference could be attributed to the internalisation of particles by cells and the difference in fluid velocity, which may increase aggregation. We also observed that, cells with internalised particles showed a different aggregation behaviour at lower flow rates while aggregation at higher flow rates follows a similar (chain) pattern as observed for ferrofluids. The aggregation behaviour observed for low flow rates might, in part, be caused by the strong field gradients at the bore of the magnet. These gradients will induce transverse motion which will increase the likelihood for aggregation as cells will come close enough to other cells to experience their magnetic dipole moment.

Cell aggregation is clearly an important factor for targeting of cells using MRT and will always be present unless cells are very dilute and or contain only a very low amount of iron oxide. The higher magnetophoretic velocity of aggregates allows for increased targeting efficiency at higher flow rates. This can be further enhanced by increasing the cell concentration in the injection solution. However, if small vessels are to be targeted, large aggregates might perturb blood flow leading to vessel occlusion. Although shear forces will increase in small capillaries, it is not certain if these forces will be big enough to disaggregate cells. However, the nanoparticles used for these experiments are super paramagnetic and the removal of the subject from the strong field of the MRI scanner would lead to the decay of dipole interactions and thus the cells should be able to disperse. This should also be the case for ferrofluids. However, there is an additional risk that during the contact time between cells, some may have bound together too strongly to be dispersed by the fluid flow. This is of particular relevance for MRT, but it will also be of importance for magnetic targeting applications that use permanent or electro magnets. Further research will be needed to address this question.

Investigation of MRT using confocal endoscopy

In order to study the effect of the imaging gradients applied to target cells we applied a gradient amplitude of 300 mT/m with a duty cycle of 2/7 perpendicular to the orientation of the fluid flow (using a gradient coil with a larger inner diameter then used for the other experiments). The supplementary movie 7 shows the effect of the gradients on the direction vector of the cells. For the case with and without gradients, 15 cells were selected at random and tracked until they left the field of view. Without any gradients applied, the direction vector for cells was 70±13°while the application of the gradients changed it to 25±14°. However we were not able to study the effect of the gradients on cell movement at lower flow rates as large aggregates rapidly blocked the field of view (cells tended to stick to the glass surface of the endoscope tip).

The influence of other parameters on targeting efficiency

Our mathematical model permits the investigation of the effects of different parameters on targeting efficiency (see section 2.3). Figure 4A for example shows the volume fraction that iron oxide particles would need to take up in a cell to reach a targeting efficiency of 90% for different gradient amplitudes.

For this plot no cell aggregation was assumed. It can be seen that such a high targeting efficiency is hard to achieve. Admittedly cell aggregation will lead to higher targeting efficiency at lower flow rates, but achieving high targeting efficiency at higher flow rates will be difficult. High targeting efficiencies will be important when MRT is used at several subsequent bifurcations to guide a bolus of cells to a target location. In order to reach the needed targeting efficiencies further optimisation will be necessary. It is clear that several factors can be influenced to maximise targeting efficiency such as: the saturation magnetisation of the superparamagnetic particles, amount of particles internalised in cells, cell concentration in bolus, viscosity, gradient amplitude and the duty cycle. The saturation magnetisation for the particles used can only be improved slightly. Different particles could be used e.g. iron cobalt nanoparticles with a high saturation magnetisation. Most of these particles are not biodegradable, limiting their potential applicability. Increasing the amount of internalised particles would be difficult as our cells are already strongly labelled and such labelling efficiencies are hard to achieve for cells with no or low phagocytic activity. The cell concentrations used in our experiments are not particularly high and can be increased easily, but the potential risk of blocking smaller vessels needs to be addressed. A lower limit for the viscosity will be given by the viscosity of an isotonic saline solution while its upper limit will be the viscosity of blood. Gradient amplitude and duty cycle can be improved via improved hardware e.g. gradient inserts optimised for targeting applications. The limit for such gradients will be cooling capacity and the maximum switching rate that does not lead to peripheral nerve stimulation.

Figure 4B shows the effect of the aggregate geometry on the targeting efficiency. As the magnetic force is a volumetric force and the drag force is dependent on the projection surface in flow direction, the deviation from a spherical object to an elongated one reduces the targeting efficiency. From the analysis of our videos we observed that the equatorial to polar diameter ratio is approximately five. If 40 cells in one spherical aggregate are necessary to explain the experimental results shown by the blue line in Figure 2A, the number of cells in one aggregate with a diameter to length ratio of 5 will increase to 60. The diameter to length ratio of aggregates is dependent on the polydispersity of the solution. It might therefore be possible to add a small concentration of ferrofluid to increase the size of these aggregates and reduce their diameter to length ratio.

Considerations for active control of cell delivery

In order to implement an active control strategy for cell delivery on a MRI scanner three steps are necessary. Firstly the current position of cells needs to be determined. Secondly the necessary gradient strength to bring them to the desired location needs to be calculated, and thirdly this gradient needs to be applied. This sequence of steps has to be repeated very rapidly to guide the cells through the vascular network as they are carried along by the blood flow. An analytical model is advantageous for the implementation of a control strategy, as calculations can be performed very rapidly. We have shown that a simple analytical model can be derived which could be used for such an application. However, the fast detection of a bolus of labelled cells in blood vessel will be challenging as the concentration of iron per millilitre of that bolus will be a few milligrams at best. The feasibility for the active guidance of a steel ball in a carotid artery of a swine has been shown previously by Martel et al [START_REF] Martel | Automatic navigation of an untethered device in the artery of a living animal using a conventional clinical magnetic resonance imaging system[END_REF]. The control strategy used relied on a pre-planned path on which the steel ball was kept using a feedback control loop which determined the difference between actual position and pre-planned position. Steering cells will be more difficult as they are not a single rigid object. An active control strategy which relies on a predictive estimate for the path cells follow could be used. This is possible with fast algebraic methods and might be useful to increase steering efficiency.

5.

Conclusions

We have presented a simple analytical model which describes the steering efficiency of magnetically labelled cells inside of a small bifurcation in the context of MRT. Such a model could be used for the implementation of an active control strategy for cell delivery using MRT. This model takes cell aggregation into account and models these aggregates as prolate spheroids. Predicted targeting efficiencies agree with experimental results when cell aggregation is taken into account, highlighting the importance of cell aggregation in this context. We have shown that changing the cell concentration or fluid viscosity in our model and experimental setup lead to the expected behaviour. We show that different parameters affecting targeting efficiency can be studied theoretically. In addition we confirmed the expected aggregation inside the MRI scanner via confocal endoscopy. Confocal endoscopy showed aggregate sizes at lower flow rates which are in the predicted range of our theoretical model. We found a transition from large aggregates containing up to 100 cells to small chain like aggregates with up to 10 cells at flow rates of 4-6 cm/s in this model. The application of imaging gradients to move cells, leads to changes of the direction vectors of cells which could be observed using near real-time confocal micro-endoscopy.

Figure 1 :

 1 Figure 1:A shows a schematic drawing of the separation area of the vascular flow phantom. An aggregate of magnetically labelled cells is symbolised by grey circles and their path over time by the fading copies of this aggregate assuming a positive X gradient. A grey circle with a dashed line indicates the field of view of the confocal endoscope. B 0 : orientation of the magnetic field in the MRI scanner, F m : magnetic force, F D : fluid drag force, R C : radius of cells, R t : radius of tube, v f : fluid velocity, a: equatorial radius of prolate spheroid as a model for aggregated cells, c: polar radius of prolate spheroid. B shows three dimensional plots for the magnetic field components (B x , B y , B z ) in the MRI scanner without a magnetic field gradient while C shows plots for the magnetic field components in the MRI scanner with a field gradient of 550mT/m applied along the x-axis. (The same coordinate system as used for the model is used for the field plots.

Figure 2 :

 2 Figure 2: A: Comparison of experimental results versus theoretical prediction for targeting efficiency. The green dots show the experimental results for a cell concentration of 4x10 6 cells/ml while the blue dots correspond to 2x10 6 cells/ml, vertical bars indicate the standard deviation. Theoretical predictions are shown as continuous lines; purple: no cell aggregation, blue: 60 cells aggregated with a prolate equatorial to polar diameter ratio of 5, green: 80 cells aggregated with the same diameter ratio. B The effect of different fluid viscosities on targeting efficiency. Blue and red dots show experimental results for PBS+3% serum and PBS+50% serum. The continuous lines show the theoretical predictions for a viscosity of 0.0016 Pa s and 0.0024 Pa s assuming no changes in aggregate size.

Figure 3 :

 3 Figure 3:Singe frame images from confocal endoscopy videos acquired at a frame rate of 12 images per second for different fluid velocities without a magnetic field and in an homogeneous magnetic field of 9.4T (inside the bore of an MRI scanner). The left column corresponds to no magnetic field while the right column corresponds to 9.4T. A: flow rate of 1 cm/s, B: 2 cm/s, C: 4 cm/s. The scale bar in the lower right hand corner of each image indicates 50 µm while the field of view is 500 µm with a depth perception of about 100 µm. The arrow in the right column figure A indicates the orientation of the magnetic field for the whole right column. D shows a schematic drawing how the microprobe was introduced into the flow phantom and how chain like aggregates are aligned with the main field of the MRI scanner (B0). In the lower right hand corner chain like aggregates of cells are magnified by a factor or three (3X).

Figure 4 :

 4 Figure 4:A Predictions for the volume fraction on iron oxide in cells necessary to reach a given targeting efficiency without cell aggregation for different gradients strengths at a duty cycle of 50%. gx: gradient strength T/m. B The effect of different polar to equatorial diameter ratio on targeting efficiency assuming the total number of cells in an aggregate is constant.

Figure 1 Figure 4 (

 14 Figure 1 (figure1a_c_new.tif)

Table 1 : Model parameters

 1

	Variable	Symbol	Value	Units
	magnetic induction	L	9.4	T
	gradient strength	ª	0.5 (0.5-10)	T m -1
	saturation magnetisation	V [	314 000	A m -1
	number of particles per cell	M	250
	tube radius	`	0.0004	m
	cell radius	`Z	0.000005	m
	average fluid velocity	b« c	0.005-0.08	m s -1
	separation length	8	0.06	m
	Duty cycle for steering gradients	Y Z	2/7
	fluid viscosity	_	0.0016-0.0024	kg m -1 s -1
	diameter length ratio	~•	5
	cells in one aggregate		1-80

Acknowledgment

The authors would like to acknowledge Carmin Habalsa (UCL) for hMNCs and Dr. Rory Warner and Dr Dan Green (Agilent Technologies, Oxford) for field maps, Prof. Roger Ordidge (UCL Biophysics) and Prof. David Hawkes (UCL CMIC) for their helpful discussions. In addition we would like to thank the British Heart Foundation, the Engineering and Physical Sciences Research Council, Guys & St Thomas' Charity, The Biological Sciences Research Council, Department of Health (DOH) UK DLH Clinician Scientist Award to Dr. Cook and GSTFT/KCL Comprehensive Biomedical Research Centre for supporting this study.

Appendix A: Cell preparation Human mononuclear cells were collected by leukapheresis from peripheral blood of G-CSF-stimulated donors. Cells were washed twice with PBS, centrifuged for 10 min at 400 g, plated in Dulbecco's Modified Eagle Medium (DMEM) (Gibco, UK) and incubated at 37 °C and 5% CO 2 for 2 hours. After two hours non-adherent cells were washed away and the media replaced with fresh DMEM + 10% fetal bovine serum (FBS) (Invitrogen, UK). Cell labelling was performed at >80 % confluency by adding 50 µL BioMag Goat anti-Rabbit IgG particles solution (Bangs Laboratories Inc, IN, USA) (particles washed 3 times with fresh media) and incubation of cells for 24 hours. After the incubation period cells were washed three times with PBS without detaching the cells. Cells were incubated with DMEM containing 25 µM CellTracker Green CMFDA (5-chloromethylfluorescein diacetate) (Invitrogen, Paisley, UK) were added and cells incubated for 30 minutes. Following that, cells were detached using PBS containing 5 mmol/l EDTA (ethylenediaminetetraacetic acid) for 10 minutes and then flushed from the surface and washed twice with PBS (by centrifugation at 300 g for 5 minutes). The labelling with the fluorescent dye was only used for experiments which included confocal endoscopy. Cells were re-suspended to 2x10 6 and 4x10 6 cell/ml in PBS containing 3% FBS or to 2x10 6 cells/ml in PBS containing 50% FBS.

Supporting Material

Videos corresponding to Figure 3: S1: no magnetic field, flow rate: 1 cm/s S2: homogeneous field, flow rate: 1 cm/s S3: no magnetic field, flow rate: 2 cm/s S4: homogeneous field, flow rate: 2 cm/s S5: no magnetic field, flow rate: 4 cm/s S6: homogeneous field, flow rate: 4 cm/s

The effect of the gradients on cell movement: S7: homogeneous field, flow rate: 3 cm/s, field gradient turned on and off Supplementary figures: S8: Supplementary figures 1 and 2