A c c e p t e d M a n u s c r i p t 1 Introduction

The topic of the presented paper is marginally related to two interesting and open problems -turbulence and the dynamo paradox. Turbulence is still considered to be one of the most challenging and not yet resolved problems. The geodynamo paradox is related to the Earth's core, in which the flows, very probably turbulent, generate and sustain the main geomagnetic field through the dynamo process. Very recently computer simulations of the geodynamo have provided models of the geomagnetic field, which give a true picture with all significant features of the observed field with its secular variations.

The paradox lies in the fact, that due to the computer limitations unrealistically high diffusion coefficients (for example viscosity) were used. Even accepting the turbulent viscosity and other diffusivities of the Earth's core, the values used in those models were noticeably higher. Thus, it was necessary to introduce hyperviscosity into the simulations to reach agreement between the observed and simulated fields. It means also that formally some type of anisotropy of diffusive coefficients was considered.

We study the conditions for marginal convection in the horizontal planar layer rotating about a vertical axis permeated by a homogeneous horizontal magnetic field. It is a model employed in [START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF], which studied both the stationary convection and the periodic instability of modes in the form of rolls inclined to the magnetic field. Considering also the case of isotropic diffusion we have advanced these studies by including two types of anisotropic diffusive coefficients. We christened the first type stratification anisotropy (SA) and it is a crude analogue of hyperdiffusivities used in computer simulations. We distinguished the diffusivities in the vertical direction from horizontal diffusivities (which are considered to be horizontally isotropic) and we introduced two subcases: oceanic type of SA and the atmospheric type of SA, respectively. In the first one (oceanic), the diffusivities in the horizontal directions are greater than in the vertical one. In the second sub-case (atmospheric) it is the reverse. The second type of anisotropy is analogous to the result obtained in [START_REF] Braginsky | Local turbulence in the Earth's core[END_REF] and it is given by the directions of the basic magnetic field, B, and the rotation, Ω . In the direction of Ω × B the diffusivities are much smaller than in the directions of Ω and B.

We focus our attention on the special type of instabilities which are affected by diffusive processes. Of course, we can expect that these diffusive instabilities, in particular, will be strongly influenced by anisotropic diffusive coefficients. Especially, if one reminds oneself that the conventional understanding of the role of diffusion -only weak damping of the instabilities which arise -is not sufficient. It is generally accepted, that the dynamics of the Earth's core are determined by three basic forces: magnetic, Archimedean and Coriolis forces (M , A, and C, respectively). However, diffusive processes may weaken the basic forces in the sense: viscosity weakens only the Coriolis force, magnetic diffusivity weakens only the magnetic force, and thermal diffusivity weakens only the Archimedean force and a new triad of basic forces, M -δM , A -δA, and C -δC, develops completely new balance of forces prone to new and even The present paper extends the papers [START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF], where the diffusive coefficients were isotropic. The paper continues in our strategy to modify old models or to set up new models of linear rotating magnetoconvection in the sense of introducing anisotropic diffusive coefficients. For instance, the results of our models of horizontal planar layer of fluid with other basic magnetic fields (homogeneous vertical or non-homogeneous azimuthal, see e.g. [START_REF] Brestenský | Linear Rotating Magnetoconvection with anisotropic diffusive coefficients[END_REF]) or with nonvertical orientation of rotation axis are upcoming for publication. We consider the horizontal magnetic field and vertical rotation axis, chosen in this paper, as a crude approximation of the Earth's azimuthal magnetic field in polar regions (not too close to poles). On the other hand, the equatorial region can be approximated by the horizontal planar fluid layer permeated by horizontal magnetic field but rotating about horizontal axis perpendicular to the field.

The results are in the form of conditions for the onset of instabilities which are set by new parameters -parameters of anisotropy of diffusive coefficients. In the case of turbulent transport phenomena the shape of transporters, i.e. eddies, distinguish transport efficiency in various directions. For that reason it is convenient to change the isotropic transport phenomena, usual in the molecular case, into the anisotropic ones. It means: for transport coefficients we have the transition from scalar into tensor quantities. Recent surveys of problems related to the Earth's core turbulence, the geodynamo paradox and some models of the geodynamo and rotating magnetoconvection with anisotropic diffusive coefficients are in [START_REF] Fearn | The geodynamo, Chapter 4[END_REF][START_REF] Schubert | Treatise on Geophysics[END_REF][START_REF] Gubbins | Encyclopedia of Geomganetism and Paleomagnetism[END_REF], where it is possible to find some important references, too.

The structure of the paper is the following: after definitions of the two types of considered anisotropies (Sect. 2), the mathematical formulation is stated and the basic dispersion relation for the two anisotropic cases is derived (Sect.

3). The method of solution for steady and overstable convection is described in Sect. 4. The influence of anisotropy on steady modes is described and numerical results are discussed in Sect. 5. Overstable modes in anisotropic cases are described in Sect. 6. The most important results are summarized and the paper is concluded in Sect. 7.

Anisotropy of diffusive coefficients

We introduce possibly the simplest forms of anisotropy, in which the thermal diffusion and viscosity are considered1 to be in the form of a diagonal Cartesian tensor with only two different values of its 3 components

ν =   ν xx 0 0 0 ν yy 0 0 0 ν zz   .

Page 3 of 27

A c c e p t e d M a n u s c r i p t

It means that in two directions the diffusion is the same, but different from the diffusion in the third direction. For instance ν xx = ν yy = ν zz or ν xx = ν yy = ν zz . We consider two basic types of anisotropy (SA and BM) of diffusive processes in the models investigated in this paper. The definitions of SA or BM type anisotropy (stratification anisotropy or Braginsky-Meytlis anisotropy) depend on the importance of basic forces determining the dynamics of the system. However, both anisotropy types SA and BM, can be formally determined by single value of diffusive coefficient, e.g. by vertical component ν zz or κ zz and by the anisotropy parameter, α.

We define the parameter of anisotropy α, a measure of anisotropy, as the ratio of two considered values of diffusivity in two different directions. Thus a measure of the anisotropy of viscosity is the parameter, α ν = ν xx /ν zz . Analogously, a measure of the anisotropy of thermal diffusivity is the parameter, α ϑ = κ xx /κ zz . Henceforth, we will assume that these parameters are equal,

α ≡ α ν = ν xx ν zz = α ϑ = κ xx κ zz , (1)
and different subscripting, if it is used, will be only for lucidity and for realising the connection between the parameters and the physical processes.

SA, Stratification anisotropy

In the SA case the density stratification determined by the vertical direction of gravity or/and the Archimedean buoyancy force has dominant influence on the dynamics of eddies. The axis of rotation has also the vertical direction z. We christen this type of anisotropy as stratification anisotropy (SA). The anisotropy is defined in the sense that diffusive coefficients have different values in the vertical direction z than in the horizontal directions x and y. We distinguish two types of this anisotropy depending on the fact that the diffusivity in the vertical direction can be greater or smaller in comparison with diffusivity in the horizontal directions (in the horizontal directions we consider isotropic diffusion, i.e. diffusion is the same in all directions of the horizontal plane, e.g. for viscosity ν xx = ν yy ). If the diffusivity in the vertical direction is greater than in the horizontal directions, we christen it a-anisotropy, which is a rough analogue of the lowermost part of the atmosphere in which turbulent transport in the vertical direction is facilitated by the developed convection. This atmospheric type of SA anisotropy is thus defined: ν zz > ν xx = ν yy and κ zz > κ xx = κ yy , and the parameter (1) of anisotropy, α, has values less than one . If the opposite case is considered, it means that the diffusivity is smaller in the vertical direction than in the horizontal directions; we call it o-anisotropy, i.e. oceanic SA anisotropy, which is analogous to the uppermost layer of the ocean, stably stratified during the day due to the heating from above. The consequence is the weakened turbulent transport in the vertical direction. In this case the anisotropy is defined by ν zz < ν xx = ν yy and κ zz < κ xx = κ yy , with the anisotropy parameter α greater than one [START_REF] Braginsky | Local turbulence in the Earth's core[END_REF].

Page 4 of 27

A c c e p t e d M a n u s c r i p t

BM, Braginsky and Meytlis type anisotropy

If rotation and the magnetic field have the dominant influence on the dynamics of turbulent eddies, we can suppose BM type anisotropy (from the study by Braginsky and Meytlis [START_REF] Braginsky | Local turbulence in the Earth's core[END_REF] of local turbulence in the Earth's core) of the system, in which the diffusive coefficients are greater in the direction of rotation (in the vertical direction) z and in the y direction of magnetic field (we suppose a homogeneous magnetic field in the horizontal direction y), than the diffusivities in the third direction x. BM type anisotropy is defined subsequently: ν zz = ν yy > ν xx and κ zz = κ yy > κ xx . The parameter (1) of anisotropy α is smaller than one in this case. In cylindrical coordinates the anisotropy results in the relationships, ν zz ∼ ν ϕϕ ≫ ν ss and κ zz ∼ κ ϕϕ ≫ κ ss , [START_REF] Braginsky | Local turbulence in the Earth's core[END_REF]. Some codes for each investigated type of anisotropy are used. Individual types of anisotropies are coded by their physical background which determine the range of anisotropy parameter values, α, too.

SA ⇒ stratification anisotropy, ν xx = ν yy = ν zz , 0 < α < 1 or α > 1
Sa ⇒ SA anisotropy of atmospheric type, a type of SA anisotropy, 0 < α < 1 So ⇒ SA anisotropy of oceanic type, o type of SA anisotropy, α > 1 BM ⇒ Braginsky and Meytlis type anisotropy, ν xx = ν yy = ν zz , 0 < α < 1

Model

Linear stability of the infinite horizontal plane layer of incompressible fluid and of the width, d, is examined. The layer is permeated by homogeneous horizontal magnetic field, B 0 , in the y direction, rotates about a vertical axis z, and gravity, g, is in the -z direction. The layer is heated from below and cooled from above. Thus the linear profile of the temperature, T 0 , is sustained; see Fig. 1.

We try to choose the parameters of the layer including diffusive coefficients so that they can correspond to the Earth's core. Therefore, we suppose that values of viscosity and thermal diffusivity are determined by turbulent eddies strongly shaped by basic forces in the core. Due to metallic electric conductivity, σ, of the Earth's core there is a good reason [START_REF] Fearn | The geodynamo, Chapter 4[END_REF][START_REF] Gubbins | Encyclopedia of Geomganetism and Paleomagnetism[END_REF][START_REF] Schubert | Treatise on Geophysics[END_REF] to neglect the turbulent contribution to the value of the magnetic diffusivity, η = 1/µ 0 σ (µ 0 is magnetic permeability). Thus η in our models is isotropic. However, our study is focused on the influence of anisotropic diffusive coefficients (coefficients of viscosity and thermal diffusivity) on the onset of thermally driven instabilities, i.e. the influence of anisotropy (1) parameter α in the SA and BM cases on rotating magnetoconvection in the form of horizontally oriented rolls.

Page 5 of 27

A c c e p t e d M a n u s c r i p t

E T © z x y -d/2 E E E E E E T B 0 E Ω = Ω 0 c T z g T 0 (z) d/2
Figure 1: Model of rotating magnetoconvection with homogeneous horizontal basic magnetic field in the infinite horizontal unstable stratified fluid layer with temperature profile, T 0 (z).

Basic equations

We consider the basic state of the system (rotating with angular velocity, Ω 0 = Ω 0 ẑ),

U 0 = 0, B 0 = B M ŷ, T 0 = T l -∆T z + d/2 d , (2)
where ŷ and ẑ are unit vectors in directions y and z, respectively. T l is a basic temperature at the bottom of the layer, z = -d/2, and ∆T is the basic temperature difference between the top and the bottom of the layer. We investigate the system stability by putting perturbations of the velocity field, magnetic field and temperature, ũ, b and θ, respectively, into the system of basic equations. This system for the unknown perturbations is converted into the dimensionless and linearized form,

R o ∂u ∂t + ẑ × u = -∇p + Λ(∇ × b) × ŷ + Rϑẑ + E z ∇ 2 α u, (3)
∂b ∂t = ∇ × (u × ŷ) + ∇ 2 b, (4)
1 q z ∂ϑ ∂t = ẑ • u + ∇ 2 α ϑ, (5)
∇ • u = 0, (6)
and

∇ • b = 0, (7)
Page 6 of 27

A c c e p t e d M a n u s c r i p t

using the following substitutions 2 , z = dz, t = d 2 η t, ũ = η d u, p = 2Ω 0 ηρ 0 p, b = B M b, θ = η∆T κ zz ϑ,
where tilded symbols are quantities with dimension and nontilded symbols are their dimensionless counterparts. In the subsequent text only dimensionless quantities are used. ρ 0 and p are density and pressure (with contribution from the centrifugal force), respectively. At first glance one can see in (3) and ( 5) that the Laplacian ∇ 2 is replaced by the anisotropic Laplacian ∇ 2 α ; see [START_REF] Gubbins | Encyclopedia of Geomganetism and Paleomagnetism[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Matsushima | The anisotropy of local turbulence in the Earth's core[END_REF]. The first equation ( 3) is the momentum Navier-Stokes equation, which describes the temporal behaviour of the perturbation velocity field u. The second equation ( 4) is the induction equation describing the temporal behaviour of the perturbation magnetic field b. The next equation of heat conduction [START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF] describes the behaviour of the temperature field perturbation ϑ in time. The last two equations express that the magnetic and velocity fields are solenoidal. In the Navier-Stokes equation we consider two forms of last term, E z ∇ 2 α u, describing the viscous force, since we consider two types of anisotropy. In the case of stratification anisotropy (SA) the viscous term has the form,

E z ∇ 2 α u = E z [(1 -α ν )∂ zz + α ν ∇ 2 ]u, (8)
and in the case of anisotropy by Braginsky and Meytlis (BM) the viscous term has the form,

E z ∇ 2 α u = E z [(α ν -1)∂ xx + ∇ 2 ]u. ( 9
)
In addition to the anisotropic coefficient of viscosity we suppose the anisotropy of thermal diffusion, too. In the equation of heat conduction we consider, therefore, two forms of the diffusive term ∇ 2 α ϑ. In SA and BM type anisotropy, respectively, the forms of diffusive terms are the following,

∇ 2 α ϑ = [(1 -α ϑ )∂ zz + α ϑ ∇ 2 ]ϑ and ∇ 2 α ϑ = [(α ϑ -1)∂ xx + ∇ 2 ]ϑ. ( 10
)
We can easily prove in both cases of anisotropy, that setting the anisotropy parameters α ν and α ϑ to one (supposing isotropic diffusivities α = α ν = α ϑ = 1), the diffusive terms in equations [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF][START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF] change into the commonly used forms, E∇ 2 u and ∇ 2 ϑ, respectively. Thus the greater α differs from one, the greater the anisotropy. We use dimensionless parameters expressing the ratio of quantities with equal dimension, for example ratios of individual forces or time scales or length scales, i.e.

R o = η 2Ω 0 d 2 , Λ = B 2 M 2Ω 0 ρ 0 µη , E z = ν zz 2Ω 0 d 2 , R = α T g∆T d 2Ω 0 κ zz , q z = κ zz η ( 11
)
2

The special form of the scaling of the temperature perturbation θ = η∆T κzz ϑ is used for the modified Rayleigh number and equation of heat conduction in their commonly used forms. By using the simpler scaling, θ = ∆T ϑ, the modified Rayleigh number has a form R = α T g∆T d/2Ω 0 η and the equation of heat conduction is then

1 qz ∂ϑ ∂t = 1 qz ẑ • u + ∇ 2 α ϑ.
Page 7 of 27

A c c e p t e d M a n u s c r i p t where R o -the modified Rossby number is a ratio of rotation period to magnetic diffusion time (the classical Rossby number R o = U/Ω 0 d is a measure of the ratio of rotation period to the time d/U , in which the characteristic flow in the system with the velocity U passes the characteristic length of the system, d); Λ -the Elsasser number is a ratio of magnetic force and Coriolis force; E z -the Ekman number is a measure of viscous forces to Coriolis force. Due to anisotropic viscosity we have also the second Ekman number,

E x = ν xx /2Ω 0 d 2 ;
R -the modified Rayleigh number measures effective buoyancy force to Coriolis force (where α T is thermal expansion coefficient); q z -the Roberts number is defined as the ratio of magnetic and thermal diffusion time scales, i.e. the ratio of coefficients of thermal diffusion and magnetic diffusion. Similarly, as in the case of the Ekman number, the Roberts number has two different values, because we consider anisotropy of thermal diffusion. Thus we have the second Roberts number,

q x = κ xx /η.
It is also useful to introduce five further dimensionless numbers

ζ = η κ zz = q -1 z , p = η ν zz = R o E z , Q = Λ E z , T a = E -2 z , R a = R E z , (12)
where ζ and p are two Prandtl numbers, dimensionless ratios of diffusivities, ζ -the inverse Roberts, and p -the inverse magnetic Prandtl number3 . The next three numbers, Q -the Chandrasekhar, T a -the Taylor, and R a -the classical Rayleigh number are useful for comparison between corresponding formulas here and in [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF][START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF].

Method of solution

In the following proceeding we split the perturbations of the velocity field u and magnetic field b into poloidal and toroidal parts,

u = a -2 [∇ × (∇ × wẑ) + ∇ × ωẑ] and b = a -2 [∇ × (∇ × bẑ) + ∇ × jẑ ]. ( 13
)
All perturbations (w, ω, b, j, and ϑ) have the form,

f (x, y, z, t) = ℜe [F (z) exp(ilx + imy) exp(λt)] , (14)
where the horizontal components of the wave vector, l and m, determine

a 2 = l 2 + m 2 .
The square root of a 2 is the horizontal wave number, a, and λ is a complex frequency (λ = iσ). The properties of convection dependent on the values l, m, and σ are discussed in the subsection 4. and ϑ) into these equations we obtain following set of ordinary differential equations for the unknown functions F (z) = W (z), Ω(z), B(z), J(z), and Θ(z),

[E z D α -R o λ] Ω + DW + imΛJ = 0, ( 15
) (D 2 -a 2 ) [E z D α -R o λ] W -DΩ + imΛ(D 2 -a 2 )B = a 2 RΘ, (16) (D 2 -a 2 -λ)J + imΩ = 0, ( 17
) (D 2 -a 2 -λ)B + imW = 0, ( 18
) (D α -ζλ)Θ + W = 0, ( 19
)
where D = d/dz and the operator D α has two different forms, that is D 2 -αl 2m 2 and D 2 -αl 2 -αm 2 in BM and SA type anisotropy, respectively. From these five equations (15-19) we successively eliminate 4 unknown functions except W and we get an ordinary differential equation of the 12th order for W of the following form,

(D α -ζλ) (D 2 -a 2 ) D λ (D α -pλ) + Qm 2 2 + T a D 2 λ D 2 + +a 2 R a D λ D λ (D α -pλ) + Qm 2 W = 0, (20)
where D λ = D 2 -a 2 -λ and the set of dimensionless parameters [START_REF] Phillips | Spherical anisotropic diffusion models for the Earth's core[END_REF] was replaced by the set [START_REF] Phillips | Spectral interactions of rapidly anisotropic turbulent viscous and thermal diffusion in the Earth's core[END_REF].

For direct comparison with [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF] and [START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF] we introduce also the substitutions,

l = l 2 π 2 , m = m 2 π 2 , σ 1 = - iλ π 2 , T 1 = T a π 4 , Q 1 = Q π 2 and R 1 = R a π 4 .
(21a) If we consider the simplest boundary conditions at z = -1/2, 1/2 (in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF] christened illustrative conditions), i.e. the stress-free and perfect thermally and electrically conducting boundaries,

W = D 2 W = DΩ = Θ = B = DJ = 0,
then due to the solution W = W 0 cos(πz) we can replace the operators in (20) by simple algebraic formulas. Finally, introducing the auxiliary variables, ã = l + m, A = 1 + ã and (21b)

A α = 1 + α l + m 1 + αã in the case of BM SA anisotropy, (21c)
we get the dispersion relation for the Rayleigh number R 1 in the form

(A α + iζσ 1 ) A [(A + iσ 1 )(A α + ipσ 1 ) + Q 1 m] 2 + T 1 (A + iσ 1 ) 2 = = R 1 (A -1)(A + iσ 1 ) {(A + iσ 1 )(A α + ipσ 1 ) + Q 1 m} . ( 22
)
The dispersion equation ( 22) for the case of the horizontal magnetic field in y direction is analogous to the vertical magnetic field case of ( [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF], chap. V).

Page 9 of 27

A c c e p t e d M a n u s c r i p t

The transition to dimensionless parameters used in [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF][START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF] serves here for comparison and checking of complex formulas, but in the following text the numerical results are presented in the space (11) of parameters R, E z , Λ, q z (or ζ = 1/q z ), p (= R o /E z ), and α.

The equation, analogous to the equation ( 22), in this more "suitable" space of dimensionless parameters has the form,

Ra 2 q z (K 2 + λ) (q z K 2 α + λ) = K 2 [E z (K 2 + λ)(K 2 α + pλ) + m 2 Λ] 2 + π 2 (K 2 + λ) 2 E z (K 2 + λ)(K 2 α + pλ) + m 2 Λ , (23)
where, using transforming formulas (21a,b,c),

K 2 = a 2 + π 2 and (24a) K 2 α = αl 2 + m 2 + π 2 αa 2 + π 2 in the case of BM SA anisotropy. ( 24b
)
In recent years, substantial progress has been made in the numerical simulations of dynamos driven by thermal convection studied here. These numerical models are governed by the same dimensionless parameters as those in the problem of rotating magnetoconvection studied here. An important limit for these models of the geodynamo is the limit of neglecting the inertial term ∂u/∂t, because it is believed that the fast modes are relatively unimportant. The basic dispersion equation in this limit of zero inverse magnetic Prandtl number, p = η/ν zz , in the space [START_REF] Phillips | Spherical anisotropic diffusion models for the Earth's core[END_REF] of dimensionless parameters E z , Λ, q z and α has the form,

Ra 2 q z (K 2 + λ) (q z K 2 α + λ) = K 2 [E z (K 2 + λ)K 2 α + m 2 Λ] 2 + π 2 (K 2 + λ) 2 E z (K 2 + λ)K 2 α + m 2 Λ . ( 25
)
Meanwhile, we stress that, in the isotropic case α = 1, there holds K 2 α = K 2 , q z = q, E z = E and similarly A α = A and D α = D 2 -a 2 , which are useful for comparing formulas (20 -25) with the corresponding formulas in [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF]. For instance, the relation (25) in the limit α = 1 is identical with (3.3) in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF].

Convection modes in the form of rolls

The main aim of linear stability analysis is to investigate the onset of the convection, thus to find the basic input parameters characterising the system and to find the configuration of instability which starts to convect most easily. The instability may occur if the Rayleigh number, the ratio of destabilising and stabilising forces, reaches a certain critical value. The task is to find the form of perturbation (determined by the wave numbers) which gives the least Rayleigh number for the basic input parameters (Ekman, Elsasser, and Prandtl numbers, and anisotropy parameter). This Rayleigh number is named the critical Rayleigh number and is labeled by R c . Searching the R c is related to the process of minimization, but first we must define all investigated modes of convection which compete for the preference in the space of dimensionless parameters [START_REF] Phillips | Spherical anisotropic diffusion models for the Earth's core[END_REF].

Page 10 of 27 A c c e p t e d M a n u s c r i p t

We investigate the onset of steady and overstable convection, which has the form of horizontal rolls and is determined by the two components, l and m, of the horizontal wave vector. The pair, l and m, can be expressed by other pair,

a = (l 2 + m 2 ) 1/2 and γ = tan -1 (m/l), ( 26
)
where a is the horizontal wave number, and γ is the angle between the axes of rolls and the imposed magnetic field in y direction. Due to [START_REF] Schubert | Treatise on Geophysics[END_REF] the vector (l, m) determines the orientation of convection rolls (rolls are perpendicular to the vector) and in the case of overstable convection this vector determines also the direction of propagating waves.

The case where the wave number l is zero we call cross-rolls, because the rolls are perpendicular to the basic magnetic field. The opposite case is when the wave number m is zero; those modes we call parallel rolls. If both wave numbers are nonzero we call these modes oblique rolls. To investigate the preference means to determine which of these forms of convection is preferred for given E z , Λ, q z , R o , and α; see [START_REF] Phillips | Spherical anisotropic diffusion models for the Earth's core[END_REF]. Combining three possible orientations of rolls with steady and overstable convection one can define 5 various modes of convection. Using codes and names by Roberts and Jones [START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF] we name them P, SC, SO, OC, and OO modes, meaning parallel, steady cross, steady oblique, overstable cross, and overstable oblique rolls, respectively. P modes are never overstable [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF].

The minimization process was done in very simple way but sufficient for this problem. Physically relevant part of the two-dimensional space of the wave vector components l and m is covered by a grid. The Rayleigh number is calculated in each point of the grid and the point with the minimal Rayleigh number is found. Suitably chosen surroundings of this point is covered by new grid with higher density of points and again the point with the minimal Rayleigh number is found. The process continues until the desired accuracy is reached with acquiring the R c and the components l c and m c of the critical horizontal wave vector. In the case of overstability formally the same minimization process was used but in each point of the grid the frequency σ was calculated. The σ is the root of the polynomial of the 7th degree which arises from the condition Im(R) = 0 where R corresponds to the expression (23). Obvious root, σ = 0, is related to the steady convection. The further roots of this equation, related to the overstable convection, stem from the polynomial equation of the 3rd degree for unknown σ 2 . In the cases of more physically possible frequencies than one, the frequency corresponding to the least Rayleigh number is critical frequency, σ c , as well as the corresponding critical l c , m c , and R c for the overstable mode. Finally, the critical Rayleigh numbers for steady and overstable modes are compared to determine the preferred mode. The determination of the orientation of the preferred modes (including also P modes with m = 0 and SC, OC modes with l = 0) is a part of the minimization process. Figure 2a shows the influence of SA type anisotropy, i.e. Sa and So on these regions. Solid and dotted lines correspond to the boundaries between individual regions in the anisotropic Sa and So cases, α = 0.5 and α = 2, respectively, and dashed lines are the boundaries in the isotropic case, α = 1. The boundary between SO and P is shifted to smaller values of the Ekman number E z in the case of Sa anisotropy. It is similar for the boundary between SO and SC, which is shifted to smaller values of the Elsasser number Λ in the Sa case. From the physical point of view we can say that Sa anisotropy enhances the effect of the Lorentz force, because the same region (region between two oblique lines -boundaries SO-SC), which in the isotropic case is named the weak magnetic field regime, corresponds in the anisotropic Sa case to the strong field regime (Λ * in the isotropic case is greater than the Λ * in the anisotropic Sa case, see the discussion in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF]). Similarly, we can say that Sa anisotropy weakens the effect of the Coriolis force, because the region between the two vertical lines, which corresponds to the fast rotation regime in the isotropic case, corresponds to a slower rotation regime in the Sa anisotropic case (E * z in the isotropic case is greater than E * z in the Sa anisotropic case). Figure 2b shows the influence of BM type anisotropy on the regime diagram. Solid lines again correspond to boundaries between the various regions for the anisotropic case (α = 0.5). The names of the regions are still valid for this BM case, even if the boundaries for the isotropic case (α = 1) are drawn for comparison by dashed lines in the figure. The boundaries between individual regions of preference are more affected by BM anisotropy than by the SA anisotropy, even if the parameter of anisotropy is the same. The change of these boundaries is not only quantitative (the boundaries are only shifted in the case of SA anisotropy) but the boundaries are changed also qualitatively because their shape is changed. The influence of BM anisotropy on the effect of the Page 12 of 27

Page 11 of 27

A c c e p t e d M a n u s c r i p t

Lorentz and Coriolis forces is thus more substantial than in the case of SA type anisotropy. Parallel modes can be preferred also in the very fast rotation regime (very small E z ) if the magnetic field is very weak (however, keep in mind that BM anisotropy is caused by stronger magnetic fields, see comment to Figure 2c). But only a small increase of the field is sufficient and its effect is considerable. The rolls which are perpendicular to the magnetic field are preferred only in a very narrow region of the Ekman number E z and the Elsasser number Λ. The distinction between fast rotation regime (E z < E * z ) and slow rotation regime (E z > E * z ) is not so simple like in the SA anisotropy case and also the isotropic case, because the boundary in the ΛE z plane between these two regimes is not simply a vertical line given by E = E * z , i.e. by the number, E * z , independently of the Elsasser number, but it is a function of it, i.e. E z = E z (Λ). It holds so that the greater magnetic field, the greater the value of the boundary Ekman number, E * z , which separates the slow rotation regime and the fast rotation regime. The third Figure 2c shows in detail the region of low Ekman and Elsasser numbers, and the region of contact of all three regimes (SC, SO, P). From the figure we can see that in the fast rotation regime arise not only SC modes but also P modes, what could be typical only for BM type anisotropy. However, it can not be the case, because BM type anisotropy must be caused by stronger magnetic fields, i.e. we must reasonably ignore the smaller Elsasser numbers, e.g. Λ ∼ < 0.01. Therefore, we focus attention only on SC and SO modes, and not on P modes in all investigated cases of anisotropic diffusive coefficients.

The critical numbers R c , a c , and γ c

Further numerical results of this stability analysis are presented in Figure 3, which express the dependencies of critical Rayleigh number R c , critical horizontal wave number a c , and critical angle of inclination of the rolls γ c on Elsasser number Λ. The main aim of this paper is to investigate the influence of the anisotropy of diffusive coefficients on the onset of rotating magnetoconvection. Therefore, various curves in the figures correspond to various values of anisotropy parameter α, which is a measure of the anisotropy. For comparison there is also the curve for the isotropic case with α = 1 in all figures.

Supposing the following range of input parameters (Elsasser number Λ ∈ (0.011, 100) and Ekman number in the z-direction, E z = 3 • 10 -7 (E x = αE z )) only two types, SC and SO, of orientation of the rolls are possible, because reasonably chosen values of E z and Λ exclude the preference of P modes (see Fig. 2). The axes of the rolls are either perpendicular to the magnetic field (angle γ c is equal 90 • ) or they are inclined to the magnetic field at some nonzero angle. These two cases are distinguished by various types of curves. For mutual comparison of the influence of the anisotropies (Fig. 3) three different values of anisotropy parameters are chosen, α = 1., 0.5, 0.1 (for completeness So is represented only by one α = 10). Various orientations of the rolls are distinguished by the type of the curves, dashed lines for cross rolls, solid lines for oblique rolls in SA anisotropy case and dotted lines for oblique rolls in the BM anisotropy case. The cross rolls are affected only by SA anisotropy, so in the BM Page 13 of 27 anisotropy case there is only one line for cross rolls, identical with the isotropic case α = 1. If we look at the equations (23, 24b) which are the formulas for the Rayleigh number R, we can see that the anisotropy parameter α is in the BM anisotropy case always multiplied by l 2 -the square of the horizontal wave number l. The horizontal wave number l is zero for rolls perpendicular to the y oriented basic magnetic field. This proves the loss of dependence of the cross rolls properties (e.g. Rayleigh number value) on BM anisotropy. Therefore, we focus our attention on the oblique rolls, i.e. on the solid and dotted lines for the anisotropy parameter values α = 0.5 and α = 0.1 (and, of course, also to isotropic α = 1).

We can see from the Figure 3a, that both BM and Sa facilitate convection by decreasing the critical Rayleigh number R c (the o type anisotropy, So with α > 1, impedes the onset of the convection, because it increases the R c ). Sa anisotropy facilitates convection more effectively than BM anisotropy, because the critical Rayleigh number is smaller in the Sa anisotropy case. This difference is the most evident for Λ = O(1) and decreases for greater values of Λ. The values of Λ for the change of preference (from cross to oblique rolls) are smaller for the BM anisotropy case. The Figure 3b shows, that the difference between critical horizontal wave numbers a c in the SA case and BM case is small. The critical horizontal wave number is slightly greater in the BM anisotropy case. This difference gradually disappears for the greater values of Λ.

We can see very clearly in the Figure 3c, that the most substantial difference between the influence of the two types of anisotropy on the magnetoconvection is the effect of anisotropies on the boundary value of Λ, for which the transition from cross to oblique rolls occurs. For example, in the case of relatively strong anisotropy (α = 0.1) the boundary values of Λ differ by one order of magnitude for the two cases of anisotropy. We can conclude, that the critical angle γ c is smaller in the BM anisotropy case and the difference between the SA and BM anisotropy cases again disappears for greater values of Λ.

Nonstationary convection, Numerical Results

From the geophysical point of view it is more interesting to analyse the influence of anisotropy on the time dependent instabilities. These instabilities are named nonstationary convection (in some papers also overstability). Nonstationary Page 14 of 27 A c c e p t e d M a n u s c r i p t convection in this model has the form of waves propagating in a certain direction (given by the basic input conditions) to the basic magnetic field (oriented in the y direction), where the angle between the magnetic field and the axes of the rolls is again one of the outputs of this study. The waves are determined by two wave numbers for two horizontal directions x and y and by frequency σ. Both magnetic field and rotation have a stabilizing effect on the fluid, because both can act as restoring forces on displaced fluid parcels and so provide a basis for possible wave motion in the fluid. Therefore, we might expect wave motions to be important in MHD systems such as the Earth's core, stars, the Sun, . . . There are many evidences for wavelike motions from measurements of magnetic field on the surface and from modeling of the field on the CMB. Waves are therefore a natural response of fluid systems to perturbations.

The Λq z regime diagrams

In the investigation of nonstationary convection new input parameters, Prandtl numbers, i.e. Roberts number, q, and inverted magnetic Prandtl number, p, appear in the system. Parameter space involves then q and p as well as Λ and E. More precisely, in anisotropy studies 4 we have p z , q z , E z and α instead of p, q, and E. The regime diagram changes from three-dimensional in the steady case to a more complicated five-dimensional in nonstationary convection. If we want to investigate the influence of anisotropy on the regime diagram of overstable modes in the simplest case, we will consider ideal electrically and thermally conducting and stress free boundaries (as an illustrative case). Further simplification is neglecting of the inertial forces and especially consideration of the often used inviscid approach (E = 0). In that case the two-dimensional Λq z regime diagram is de facto a cross section of the general three-dimensional diagram ΛE z q z in the plane E z = 0. In the following Figure 4 there are the two Λq z regime diagrams derived numerically. Some important details in the diagrams have been confirmed asymptotically, but we do not present it.

The Figures 4ab show the influence of two anisotropies on this diagram. The curves on these figures represent boundaries between areas with a preferred type of orientation of the rolls. We consider like in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF] The influence of the stratification anisotropy is described in the Figure 4a. Various types of curves correspond to different cases of SA anisotropy (different values of the anisotropy parameter α) and the isotropic case has solid curves. The boundary between the SO and OO regimes is vertical line independent of the Elsasser number Λ and it is in essence a concrete number, q z = q * z . This boundary is not an existence boundary of overstable modes, but the boundary of preference. The boundary between the SO and SC regimes is a line, Λ = Λ * (α), independent of q z (steady convection is always independent of the Roberts number), where Λ * (α) is an increasing function of α. This fact was discussed with Fig. 2 for E z = 0 in the case of the ΛE z regime diagram for steady convection, too. Boundary OO-OC is again a line representing linear dependence Λ = Λ(q z ) with inclination increasing with α. The remaining boundaries between regimes SO-OC and SC-OC are complicated curves, Λ = Λ(q z ), and their analytical derivation is also complicated. The main task is to investigate how these boundaries are changed under the influence of the anisotropy of the diffusive coefficients. At first sight one can see that SA anisotropy changes the boundaries only quantitatively -the basic nature is preserved. The whole set of boundaries for fixed α is shifted towards the smaller Roberts numbers and greater Λ in the oceanic type anisotropy (So, α > 1). In the case of atmospheric type (Sa, α < 1) it is shifted towards the greater Roberts numbers and smaller Λ.

The second Figure 4b describes the influence of BM type anisotropy on the regime diagram. Solid lines are related to the isotropic case and the remaining anisotropic cases are distinguished by different types of curves. Comparing with the previous Figure 4a one can see that BM anisotropy changes the regime diagram not only quantitatively but even qualitatively. The main difference is related to the SO-OO boundary, which is a simple vertical line q z = q * z independent of the Elsasser number only in the isotropic case. Even weak anisotropy (α = 0.9) changes it into the curve Λ = Λ(q z ). The obliqueness increases with increasing anisotropy. The next important difference from the SA case is connected with the OC-OO boundary, which is independent of the anisotropy parameter, the same like the isotropic case. The boundary between SO and SC is similarly a horizontal line like from the SA anisotropy case in the previous Figure 4a. However, this similarity of boundaries SO-SC in SA and BM type anisotropies is not obvious in non-zero E z cases as is seen in Fig. 2, showing the regime diagram in the ΛE z -plane for steady convection. The most probable reason is in the well known big difference between the case, E → 0, and the here considered case, E = 0. The remaining boundaries between the regimes SO-OC and SC-OC are again complicated curves Λ = Λ(q z ) as in SA anisotropy. The Figures 5ab show γ c and σ c as functions of Λ for the high rotation case, E z = 3 • 10 -7 , in BM type anisotropy (R c and σ c are not presented in figures, see comment below.). Seven curves refer to seven various values of the anisotropy parameter α = 1, 0.7, 0.6, 0.55, 0.5, 0.3, 0.1. Obviously, at α values, 0 < α < 1, typical for BM anisotropy, it is valid that the smaller anisotropy parameter α, the stronger the anisotropy, i.e. at α → 0 the BM anisotropy is the strongest one, and at α → 1 is the weakest one. At first sight a considerable effect of the anisotropy is seen in particular in the reduction of the existence interval of Λ in which the instabilities can occur. The interval decreases from the right side and the left side is independent of anisotropy. This means that even strong BM anisotropy can not change the fact that overstable modes exist only for sufficiently strong magnetic field (i.e. the Elsasser number has to be large enough; see the left-hand end of the curves). In the first Figure 5a one can see that nonstationary convection starts to exist as cross rolls (γ c = 90 • ) at a certain value of Elsasser number, but for greater values of Λ the rolls are more oblique to the direction of magnetic field. The frequency is independent of the Elsasser number only in the isotropic case, but in anisotropic cases the frequency is a decreasing function of Λ with a sharp drop to zero, which is seen in the Figure 5b. The value of Λ for this sharp drop is an increasing function of the anisotropy parameter α, meaning that this value is smaller for the stronger BM type anisotropy. This limitation of overstable modes existence is also evident from the Figure 5a as the sharp right-hand end of the curves. This is one of the important discovered influences of BM type anisotropy, because there is no indication for this property in the isotropic case or in SA type anisotropy. The critical Rayleigh number R c and critical horizontal wave number a c as functions of Λ are not shown, because there is no influence on these dependencies in the Page 17 of 27 A c c e p t e d M a n u s c r i p t case of BM anisotropy for very small E z = 3 • 10 -7 as well as in the isotropic case [START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF]. R c is not dependent on Λ (and nor α for OO -overstable oblique modes), because R c = 6ζ √ 3π 2 ( . = 51.28 for ζ = 0.5 where ζ is q -1 z ). Similarly, the critical horizontal wave number a c has the same feature. Moreover, a c = π √ 2 . = 4.44 is constant and even it does not depend on ζ like R c . There is another reason to not show the dependencies of R c (Λ) and a c (Λ). Varying length curves of these dependencies at different α for BM anisotropy overlap and therefore are indistinguishable. Next Figure 6 shows the same dependencies (now also including R c and a c as functions of Λ) again in the fast rotation case, but now for higher values of the Ekman number E z = 10 -3 , but for the same Prandtl numbers p = 0.9 and ζ = 0.5 like in Fig. 5. Various curves correspond to values of the anisotropy parameter α = 1, 0.7, 0.6, 0.55, 0.5, and 0.3, as for smaller E z = 3 • 10 -7 in this BM case.

Page 16 of 27

The effect of BM anisotropy is obvious from all four Figures 6abcd and the differences between isotropic and anisotropic cases increase more strongly with increasing anisotropy than for the much smaller Ekman number E z = 3 • 10 -7 case. Equally as in the preceding case of smaller E z (Fig. 5) BM anisotropy affects non favourably the onset of nonstationary convection from a particular point of view, because it narrows the existence interval of the Elsasser number. From the Figure 6a one can see that the minimal value of Λ for existence of overstability slightly decreases, but the critical Rayleigh number increases for increasing anisotropy. The critical horizontal wave number a c has similar behaviour; it is a slightly increasing function of BM anisotropy. Critical frequencies are decreasing functions of the Elsasser number and have small values. The sharp drop of frequencies and subsequent end of the existence of overstability for higher values of the Elsasser number is seen for all the considered values of the anisotropy parameter α. The value of Λ at which the sharp drop of frequency occurs is a decreasing function of BM anisotropy (the smaller anisotropy parameter α, the smaller Λ).

The same dependencies were investigated also in the case of stratification anisotropy. The influence of this SA type anisotropy on nonstationary convection is very simple and therefore is not presented in figures. The critical Rayleigh number R c as well as the critical horizontal wave number a c are not influenced by SA type anisotropy and have the same value as in the isotropic case [START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF] The next Fig. 7 is again related to the BM case, but shows R c , a c , σ c and γ c as functions of the Ekman number E z for the constant Elsasser number Λ = 3. The Prandtl numbers p = 0.9 and ζ = 0.5 are the same like in other nonstationary studies presented in Figures 5 and6.

The critical Rayleigh number R c is for smaller Ekman number E z not dependent on it and is not influenced by BM anisotropy. The slight dependence of R c on E z and α is only for the greater E z > O(10 -5 ). For smaller Ekman numbers the critical horizontal wave number is not influenced either by the anisotropy or by the Ekman number alone. However, for Ekman numbers greater than E z = O(10 -5 ) the influence of anisotropy starts to be evident and the horizontal wave number behaves as an increasing function of the Ekman number for the case of strong anisotropy (α = 0.01) as well as a decreasing function of the Ekman number for weak anisotropy up to the isotropic case (α = 1), with a smooth transition between these two behaviours. We can see similar behaviour in the Figure 7d where the critical angle γ c is a function of anisotropy but independent of the Ekman number for small Ekman numbers. For greater E z the dependence on E z itself appears and the effect of BM anisotropy is also qualitative. In the isotropic case the angle γ c increases with increasing Ekman number so the rolls tilt into the direction perpendicular to the magnetic field. BM anisotropy changes this behaviour into the opposite one -for greater Ekman numbers the rolls tilt into the direction parallel to the magnetic field (the angle γ c decreases). The third Figure 7c shows that sharp drop of frequency occurs for the cases of stronger anisotropy, the nonstationary convection does not exist for values of the Ekman number greater than the value at which the frequency reaches zero (this value is an increasing function of the anisotropy parameter α). The BM type anisotropy handicaps the overstability by reducing the range of Ekman numbers at which the overstability exists (similar to the case of the Elsasser number).

Page 19 of 27 A c c e p t e d M a n u s c r i p t

Conclusions

The influence of anisotropy of the diffusion coefficients (thermal diffusion and viscosity) on the marginal stability of a horizontal planar layer of electrically conducting fluid has been studied. The layer is permeated by a horizontal homogeneous magnetic field and rotates about a vertical axis of rotation. This model was employed by Roberts and Jones [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF] and also by Eltayeb [START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF]. Considering yet the isotropic diffusion we advanced these studies by including two types of anisotropic diffusive coefficients, SA -stratification anisotropy (ν zz = ν xx = ν yy , κ zz = κ xx = κ yy ) and BM -the anisotropy of Braginsky and Meytlis [START_REF] Braginsky | Local turbulence in the Earth's core[END_REF]

(ν zz = ν yy > ν xx , κ zz = κ yy > κ xx ).
Three regimes, as in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF], are possible in the case of steady convection: P -rolls parallel to the magnetic field; SO -oblique rolls; and SC -cross rolls. Both types of anisotropy influence the ΛE z regime diagram (Fig. 2). The boundaries between individual regions of preference are more affected by BM anisotropy than by SA anisotropy in comparison with the isotropic case. SA anisotropy of atmospheric type enhances the effect of the Lorentz force and weakens the effect of the Coriolis force. In BM anisotropy the change of these boundaries is not only quantitative (the boundaries are only shifted in the SA anisotropy case), but the boundaries are changed also qualitatively because their shape is changed. The influence of BM anisotropy on the effect of the Lorentz and Coriolis forces is thus more substantial than in the SA anisotropy case. SA anisotropy of atmospheric type facilitates convection more effectively than BM anisotropy, because the critical Rayleigh number is smaller in this SA case. Critical horizontal wave numbers are slightly greater in the BM anisotropy case than in the SA and isotropic cases. The differences between the influence of SA and BM anisotropy are the most evident for Λ = O(1) and decrease for greater values of Λ.

Supposing nonstationary convection new input parameters appear. Neglecting inertial forces and considering an inviscid approach the influence of both types of anisotropy on the Λq z regime diagram has been studied. Four regimes are compared, as in [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF][START_REF] Eltayeb | Hydromagnetic convection in a rapidly rotating fluid layer[END_REF][START_REF] Eltayeb | Overstable hydromagnetic convection in a rapidly rotating fluid layer[END_REF]: steady cross rolls (SC); steady oblique rolls (SO); overstable cross rolls (OC), and overstable oblique rolls (OO). SA anisotropy changes the boundaries between regimes only quantitatively -the basic nature is preserved (in comparison with the isotropic case [START_REF] Roberts | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer I. Finite Magnetic Diffusion[END_REF]). The region boundaries are shifted towards the smaller Roberts numbers and greater Elsasser numbers in the case of oceanic type anisotropy (α > 1), and in atmospheric type (α < 1) they are shifted in opposite directions. BM anisotropy changes the regime diagram not only quantitatively but even qualitatively. In BM anisotropy the OC-OO boundary is independent of the anisotropy parameter and is the same as for the isotropic case.

The influence of anisotropies on critical numbers (R c , a c , γ c and σ c ) as functions of the Elsasser number Λ and the Ekman number E z has been also investigated. In the case of very small E z = 3 • 10 -7 the critical Rayleigh number R c is not dependent on Λ (and nor α for OO -overstable oblique modes), because R c = 6ζ The effect of BM anisotropy is more substantial for larger Ekman number E z = 10 -3 (Fig. 6). The difference between the isotropic and anisotropic BM cases increases more strongly with increasing anisotropy than for the much smaller Ekman number E z = 3 • 10 -7 case. R c and a c are no longer simple constants, but both are functions of Λ and the anisotropy parameter α.

There is one unique property of BM anisotropy (in comparison with the SA or/and isotropic cases). BM anisotropy handicaps the overstability by reducing the ranges of the Elsasser number as well as the Ekman number at which overstability occurs (see e.g. drops of frequencies in Figs 5 -7 in curves σ c (Λ) and σ c (E z )). This reduction is in the sense that greater the BM anisotropy, the smaller the maximum Λ and E z at which nonstationary convection ceases to occur.

We expect that additional study (in preparation) of anisotropic diffusivities influence on rotating magnetoconvection matched to equatorial regions will reveal other interesting facts.

Page 21 of 27

1 .

 1 Taking the z-components of the curl and double curl of the Navier-Stokes equation, the curl of the induction equation and the induction equation itself and the equation of heat conduction and substituting the perturbations (w, ω, b, j Page 8 of 27 A c c e p t e d M a n u s c r i p t

A c c e p t e d M a n u s c r i p t 5

 5 Stationary convection, Numerical results 5.1 The ΛE z regime diagrams The following numerical results (Fig.2) show the influence of the anisotropy on the distribution of the regions of preference in the ΛE z -plane, where three regions are coded by the orientation of the preferred rolls P ⇒ Steady Parallel roll is preferred; SC ⇒ Steady Cross roll is preferred; SO ⇒ Steady Oblique roll is preferred.

Figure 2 :

 2 Figure 2: The regime diagram for steady convection in cases of both type of anisotropy, SA and BM. Regime diagram shows the regions in the ΛE z -plane with the preferred type of convection.

Figure 3 :

 3 Figure 3: Comparison of both types of anisotropy. Dependence of a) critical Rayleigh number R c , b) critical horizontal wave number a c and c) critical angle between rolls and basic magnetic field γ c , on Elsasser number in the case of stratification anisotropy (SA) and anisotropy by Braginsky and Meytlis (BM). Cross rolls -dashed curves, oblique rolls -solid curves (SA), dotted curves (BM).

Figure 4 :

 4 Figure 4: Regime diagrams in the cases of both types of anisotropy, a) SA and b) BM. Regime diagrams show the regions in the Λq z plane with the preferred type of convection in the inviscid case, E z = 0. (The position of symbols SO, SC, OO and OC corresponds to correct regions bounded by solid lines for isotropic case.)

  four regimes SC ⇒ Steady Cross roll is preferred; Page 15 of 27 A c c e p t e d M a n u s c r i p t SO ⇒ Steady Oblique roll is preferred; OC ⇒ Overstable Cross roll is preferred; OO ⇒ Overstable Oblique roll is preferred.

A c c e p t e d M a n u s c r i p t 6 . 2

 62 The dependence of R c , a c , σ c , and γ c on Λ Similarly to the case of steady convection our attention is focused on the influence of anisotropy on the dependencies of the basic parameters of instability (critical Rayleigh number R c , critical horizontal wave number a c , critical angle γ c and critical frequency σ c ) on input parameters decribing the basic state (Elsasser number Λ, Ekman number E z and anisotropy parameter α, but for fixed values of the Prandtl numbers, ζ and p). The following results are presented in the form of figures of dependencies of R c , a c , σ c and γ c on the Elsasser number for constant Ekman number (or dependencies on the Ekman number with constant Elsasser number), when the inertial forces are considered (non-zero Rossby number). Various strengths of anisotropy are represented by various curves in the figures (curves are labeled by value of the anisotropy parameter α).

Figure 5 :

 5 Figure 5: The influence of BM type anisotropy on nonstationary convection. a) critical angle between rolls and basic magnetic field γ c and b) critical frequency σ c as functions of Elsasser number Λ for constant Ekman number E z = 3 • 10 -7 .

Figure 6 :

 6 Figure 6: The influence of BM type anisotropy on nonstationary convection. a) Critical Rayleigh number R c , b) critical horizontal wave number a c , c) critical frequency σ c and d) critical angle between rolls and basic magnetic field γ c as functions of Elsasser number Λ for constant Ekman number E z = 10 -3 .

(R c = 6ζ √ 3π 2 where ζ is q - 1 z and a c = π √ 2 . = 4 .Figure 7 :

 1247 Figure 7: The influence of BM type anisotropy on nonstationary convection. a) Critical Rayleigh number R c , b) critical horizontal wave number a c , c) critical frequency σ c and d) angle between rolls and basic magnetic field γ c as functions of Ekman number E z for constant Elsasser number Λ = 3.

√ 3π 2 2 . = 4 .

 224 ( . = 51.28 for ζ = 0.5). Similarly, the critical horizontal wave num-Page 20 of 27 A c c e p t e d M a n u s c r i p t ber a c has the same feature. Moreover, a c = π √ 44 is constant. These features hold for both anisotropies. It means, that anisotropy does not change the preference of the overstable modes in the very small E z = 3 • 10 -7 case.

1 α

 1 Figure(s)

see expressions[START_REF] Gubbins | Encyclopedia of Geomganetism and Paleomagnetism[END_REF][START_REF] Jones | The Onset of Magnetoconvection at Large Prandtl Number in a Rotating Layer II. Small Magnetic Diffusion[END_REF] and[START_REF] Matsushima | The anisotropy of local turbulence in the Earth's core[END_REF] for anisotropic Laplacian ∇

α in equations (3) and (5), respectively.

Due to the next ignoring of pressure we use symbol p for the one frequently used dimensionless number, p = pz.

In subsection 3.1 we have already defined[START_REF] Phillips | Spectral interactions of rapidly anisotropic turbulent viscous and thermal diffusion in the Earth's core[END_REF] anisotropic Prandtl numbers p = pz and ζ = q -1 z .

Acknowledgements

This work was supported by scientific grant agency VEGA (grant n. 1/3069/06) and by Comenius University (grant n. UK /398/2008). We are also grateful to Sebastián Ševčík and Alexandra Marsenić for helpful discussions and hints. We would like to thank also to the organizers of SEDI meeting in Kunming for financial support to present our results by TS. We are also happy to express our gratitude to both unknown reviewers for very important comments and helpful hints, which has helped us to revise our manuscript to be more clear and readable.


C * ; Madeira

Sc ; Ribeiro

M Z ; Turk

J.M.S. Cabral

Optimization of gene delivery to HEK293T cells by microporation using a central composite design methodology

Keywords: Cell recovery, Central Composite Design, Gene delivery, Microporation, Plasmid DNA

Microporation is an efficient method for delivering plasmid DNA molecules into cultured cells.

Herein, we present the optimization of gene delivery by microporation using a Central Composite Design methodology. It was given relevance not only to the transfection efficiency but also to the cell recovery. Different amounts of DNA (1 and 3 g) mainly affected cell viabilities and cell recoveries, which decrease from 93 % to 76 % and from 47 % to 25 % respectively, when higher DNA quantity is used. With this work we suggest an easy methodology to improve transfection of mammalian cells underlining the feasibility to achieve 60% of gene delivery efficiencies whilst recovering 50 % of cells, with 90 % of viability.

Introduction

Gene electrotransfer has been a powerful tool to transfer DNA molecules into a large variety of prokaryotic and eukaryotic cells. However, there are still some problems urging to be solved such as low cell viabilities and low transfection efficiencies especially when using low DNA amounts.

Most of published data has been obtained from conventional electroporation that usually relies on the use of high electrode surfaces areas requiring minimum volumes of 50 l. Novel devices with different layouts have been recently developed and tested with several cell lines for diverse applications (see recent review [START_REF] Lee | Microscale electroporation: challenges and perspectives for clinical applications[END_REF])). Our study comprises the use of a microcapillary system with an electrode surface area of 0.33 mm 2 and 10 l chamber for cell suspension microporation [START_REF] Kim | A novel electroporation method using a capillary and wire-type electrode[END_REF]. In most published work, transfection efficiency and cell viabilities are the only considered output values for electroporation optimization, which we consider that may mislead the process efficiency. In fact, high cell viabilities may be obtained even when few cells are recovered after being submitted to electric voltages. To our best knowledge few authors have published results regarding cell recovery and yield of transfection after electroporation of mammalian cells [START_REF] Aluigi | Nucleofection is an efficient nonviral transfection techni que for human bone marrow-derived mesenchymal stem cells[END_REF][START_REF] Ferreira | Optimization of a gene electrotransfer method for mesenchymal stem cell transfection[END_REF].

In this work a Central Composite Design (CCD) methodology was used as a valuable tool to assess the effect of several independent factors on HEK293T gene delivery efficiency by microporation. These cells are able to replicate appropriated plasmids [START_REF] Wurm | Transient gene expression: A novel mammalian cell-based technology for recombinant protein production -A historical and technical perspective[END_REF] and have been used in basic research and for production of virus and recombinant proteins [START_REF] Pham | Large-scale Transfection of mammalian cells for the fast production of recombinant protein[END_REF].

Moreover, microporation is a recent available and promising technology of gene transfer and taking into account the conditions that gave rise to higher percentages of cells expressing a reporter protein we have evaluated not only cellular viability but also cell recovery and transfection yield which constitutes new data regarding this methodology. Taken together, our results will certainly be very useful in gene delivery applications to mammalian cells.

Materials and methods

Plasmid construction and purification

The plasmid pcDNA3.1 (Invitrogen, La Jolla, CA) encoding green fluorescent protein (GFP) reporter gene was constructed using standard cloning methods. This plasmid has a SV40 ori, therefore able to replicate in mammalian cell lines that stably and constitutively expresses the SV40 large T antigen such as HEK293T. In brief, the reporter gene was digested from the plasmid pEGFP (BD Biosciences Clontech, Mountain View, CA) using restriction enzymes Hind III and Apa I and further cloned into pcDNA3.1. The construct (5,7 kb) was produced in E. coli and further purified using the Endotoxin-free (EF) Plasmid DNA Purification, NucleoBond ® Xtra Midi EF (Macherey-Nagel, Düren, Germany) and resuspended in EF-water at 1 g/ l.

Cell culture and microporation Human Embryonic Kidney (HEK) 293T cells (ATCC-LGC Nr: CRL-11268; Middlesex, UK) were grown in adherence in T 75 flasks with Dulbecco's Modified Eagle Medium (DMEM) high glucose (4.5 g/l), 10% fetal bovine serum (FBS), 1% penicillin/streptomycin (100 μg/ml) (Gibco, Invitrogen, La Jolla, CA) and were maintained at 37 °C and 5% CO 2 -humidified atmosphere.

Sub-confluent cells were harvested and washed with DPBS. Cells were resuspended in a sacarose based buffer (SB-PBS: 250 mM Sacarose and 1 mM MgCl 2 in DPBS) or in the Resuspension Buffer (RBprovided by the equipment manufacturer) at a density of ~2 × 10 5 cells/10 l and incubated with plasmid DNA on ice. Microporation was accomplished within 20 min. maximum using a Microporator (MP-100) purchased from Digital Bio Technology (Seoul, South Korea) according to manufacturer's instructions. The ratio between cells and DNA ranged 7 10 4 -20 10 4 cells/ g and DNA amount is referred when appropriated. Afterwards cells were added to 100 l of pre-warmed medium without antibiotics and equally divided into 2 wells of a 24-well culture plates containing culture medium under the same conditions. Cells were incubated at 37 °C in a 5% CO 2 -humidified atmosphere, harvested after 24 hours using DPBS and analyzed.

Transfection efficiency, Cell Viability, Cell Recovery and Yield Transfection efficiency was measured by flow cytometry using FACSCalibur equipment (BD Biosciences, San Jose, CA) and reported as percentage of GFP positive cells (GFP + cells). A laser excitation wavelength of 488 nm was used and experiments were conducted in 137 mM NaCl, 2.8 mM KCl solution as the sheath fluid. Results considered statistically representative had minimum 1000 events gated. Data were analyzed by CellQuest software (BD Biosciences, San Jose, CA).

The percentage of viable cells was estimated by trypan blue exclusion method. For each micro-electroporated sample (m), Cell Recovery (CR) was determined using the equation:

100 (%) c m m CA CA CR (Eq.1)
where CA is the number of cells alive and c is the non-electroporated control. The yield of transfection (Y) was calculated using the following equation:

c m m CT CA GFP Y (%) (Eq.2)
where CT is the number of total cells.

Experimental Design

Factorial designs are generally used to evaluate the effects of multiple independent factors and their degree of interaction. In this work we used a CCD to optimize gene delivery to HEK293T cells, using a microporation technology. The response (dependent variable) was the percentage of GFP + cells measured by flow cytometry. Accordingly, we evaluated the effects and interaction of four independent variables: electric pulse (X 1 ), pulse width (X 2 ), number of pulses (X 3 ) and DNA amount (X 4 ). Using a face centred star points design, three coded levels (low (-1); central (0) and high (+1)) were chosen according to previously obtained results (data not shown). The coded values of independent variables are given in Table 1.

(Table1)

This CCD design was used with the assistance of STATISTICA software (StatSoft, Tulsa, OK, USA) which has allowed estimation of a full quadratic model with the general description: N = 2 k-p + 2k +C 0 , where N is the number of experiments, k is the number of independent variables (k = 4), p the fractionalization number (in a full design, p = 0) and C 0 is the number of central points (C 0 = 4), required for curvature estimation [START_REF] Barker | Statistical modelling of the formulation variables in non-viral gene delivery systems[END_REF]. Accordingly, a total of 28 (2 4-0 +2 4+4) different combinations (including four replicates of the central point each sighed the coded value 0) were chosen according to a CCD configuration for four independent variables.

Results and Discussion

The majority of reported studies regarding optimization of gene delivery efficiency to [START_REF] Barker | Statistical modelling of the formulation variables in non-viral gene delivery systems[END_REF]. Similarly other authors have verified the effect of several parameters and their interaction on physical stability of polyplexes [START_REF] Gazori | Evaluation of Alginate/Chitosan nanoparticles as antisense delivery vector: Formulation, optimization and in vitro characterization[END_REF][START_REF] Zhong | Preparation and characterization of a novel nonviral gene transfer system: Procationic-liposome-protamine-DNA complexes[END_REF]. Regarding electroporation, and to our knowledge there is no reported study aiming at multivariate optimizing procedure for in vitro gene delivery in mammalian cells.

With our work we intended to show how the variables: electric pulse (X 1 ), pulse width (X 2 ), number of pulses (X 3 ) and DNA amount (X 4 ) affect gene delivery by microporation, and also inquire about their main interactions using a Central composite design (CCD). The number of GFP + cells, obtained by flow cytometry analysis was the dependent variable and the results are shown on Table 2. When using an electric pulse of 1500V within 30ms and/or 3 pulses, the recovered cells were not enough for flow cytometry measurement due to high cell mortality and so we considered as zero in the CCD.

(Table 2)

A second-order polynomial equation (Eq. 3) was used to express the GFP + cells percentage as a function of four independent variables. The quadratic effect of Electric pulse (X 1 2

) and DNA amount (X 4 2

) have shown the most pronounced effects on the number of GFP + cells.

GFP + cells (%) = 36.1+ 2.8X 1 -3.1X 2 -1.8X 3 + 2.7X 4  25.7X 1 2 -9.0X 2 2  4.4X 3 2 + 9.6X 4 2  6.3X 1 X 2 -3.6X 1 X 3 + 1.6X 1 X 4 + 3.4X 2 X 3  1.7X 2 X 4 (R 2 =0.95) (Eq. 3)
For each parameter estimator a test of significance was performed. Parameters with less than 95% significance (p-values above 0.05) were discarded and pooled into the error term (residual error). Accordingly, in the previous equation the parameter linear relation of X 3 by X 4 was not considered (p = 0.56). The model was further validated by analysis of variance (ANOVA), using a methodology previously presented by others [START_REF] Martins | Application of central composite design for DNA hybridization onto magnetic microparticles[END_REF], for determining which of the factors significantly affect the response variables in study (GFP + cells), using a Fisher's statistical test (F-test),. The results are shown on Table 3. (Table 3)

Data show that the model is statistically significant (for all the regression models, F-value > F-critical, p-value < 0.001) explaining 95 % (R 2 = 0.95) of the observed variance. There is no evidence of lack of fit at 95 % confidence level, as determined by the ratio MS(LF)/MS(PE) (for all response variables, F-value (LF) < F-critical and p-value > 0.05), which means that the model explains the observed differences on the response variable. The relationships between independent and dependent variable can be easily visualized using RSM, where the obtained regression model was used to calculate the response surface for each response variable (Fig. 1).

(Figure 1)

The effects of electric pulse and DNA amount on GFP + cells percentage are visualized from a threedimensional plot and a contour plot (Fig. 1A). A saddle-shaped surface plot was obtained and two maximal points were observed: one at higher and other at lower DNA amount, using an electric pulse of ~1200 V. Conversely, when using two pulses and 2 g of DNA the highest percentage of GFP + cells are obtained at electric pulse and pulse width ranging from 1200-1300 V and 20-24 ms respectively (Fig. 1B). A similar peak-shaped response surface was obtained when plotting electric pulse and number of pulses at 20 ms and using 2 g of DNA, and again a higher number of GFP + cells was obtained with 1200 V (Fig. 1C). Optimal experimental conditions were selected from the response surfaces plots and higher values of GFP + cells (~60 %) were obtained with 1200-1300 V, 20 ms of pulse width either using 1 or 3 g of DNA. The number of pulses (1-3) did not greatly affect the number of GFP + cells contrarily to some reported studies, although using different cell lines and electroporation equipments [START_REF] Lin | Efficient expression of foreign genes in CHO DHFR cells by electroporation[END_REF]. It is well known that when cells are submitted to high electric pulses a major amount of DNA enter into the cells though concurrently an increase of cell death is observed [START_REF] Ferreira | Optimization of a gene electrotransfer method for mesenchymal stem cell transfection[END_REF]. Several studies, using different cells, have shown that a way to circumvent this fact is by decreasing electric pulse and increasing the duration of pulse (pulse width) or within low voltages use more than one pulse [START_REF] Cemazar | Control by pulse parameters of DNA electrotransfer into solid tumors in mice[END_REF]). Higher transfection efficiency levels have also been obtained by increasing the DNA amount up to 60 g (corresponding to ~2 10 4 cells/ g DNA) [START_REF] Helledie | A simple and reliable electroporation method for human bone marrow mesenchymal stem cells[END_REF]) regardless lower DNA amounts are commonly employed. In fact, in our experiments with 3 g of DNA a similar ratio of cell/ g DNA (~7 10 4 ) was used.

We further compared two different microporation buffers using different DNA amounts with the purpose of evaluating respective cellular viabilities, cell recoveries and overall yield. In this study we took advantage from the information retrieved from the previously shown CCD experiment and fixed the following microporation conditions: 1250 V, 20 ms and 1 pulse. One buffer was supplied by the equipment kit (RB -with unknown formulation) and the other was prepared in the laboratory with the formulation: 250 mM Sacarose, 1 mM MgCl 2 and 10 mM Tris.HCl in DPBS (SB-PBS). The data shown in Figure 2 were obtained from three independent experiments. (Figure 2)

Using the equation model obtained by the experimental design (Eq.3) the predicted values for GFP+ cells when using SB-PBS under the referred conditions (in coded level: X 1 = 0.165; X 2 = 0; X 3 =  1) are 40 % for 1 g DNA (X 4 =  1) and 46 % for 3 g DNA (X 4 = 1). In fact, from figure 2A the obtained values were respectively 55% and 63%, slightly higher than the predicted values, though as predicted both had similar values. On the other hand, when using RB the DNA amount is crucial. Lower amount of DNA displays higher number of GFP + cells (Fig 2B). As depicted from Fig. 2 the buffer prepared in the laboratory (SB-PBS) displayed higher cellular viabilities, and most importantly higher cellular recoveries and higher yields. Common to both buffers, as expected, is the fact that lower DNA amounts give rise to higher cell viabilities and recoveries which is mainly reflected in higher yield. Cellular viability using SB-PBS ranged from 70-90 % while with RB it achieved a maximum of 60 % using a lower DNA amount. Cell recoveries and transfection yields obtained with RB were lower than 20 % whereas with SB-PBS a maximum recovery of 50 % and 30 % yield were obtained. These results suggest that RB of unknown composition might have a higher conductivity than SB-PBS. In fact, previous reports have shown that the percentage of surviving cells increases with the decreasing medium conductivity [START_REF] Ferreira | Optimization of a gene electrotransfer method for mesenchymal stem cell transfection[END_REF] suggesting that electroporation buffer is other parameter that might be considered on electrotransfer optimization protocols. Even though each electroporation equipment has a specific buffer suggested by the supplier, some authors have recently reported promising results of gene delivery efficiencies using accessible buffers of known composition [START_REF] Ferreira | Optimization of a gene electrotransfer method for mesenchymal stem cell transfection[END_REF][START_REF] Kang | Phosphate-buffered saline-based nucleofection of primary endothelial cells[END_REF]. Actually, in one experiment transfection efficiency obtained with RB was slightly higher than those obtained with SB-PBS because 78 % of GFP + cells were obtained using 1 g of DNA. Contrarily, increasing DNA amount led to a clear reduction of transfection efficiency. This fact is in fair agreement with manufacturer instructions which suggests the use of 0.5 g of DNA per 1.5 10 5 cells. It was recently reported that electroporation of cells in a capillary system gave rise to an enhancement of transfection rate and cell viability when compared to the cuvette system [START_REF] Kim | A novel electroporation method using a capillary and wire-type electrode[END_REF]). In the conventional electroporation, using a cuvette, the reported cellular viabilities are generally below 70 % and most importantly the cell recoveries are rarely mentioned, with the exception of few studies with primary cells [START_REF] Aluigi | Nucleofection is an efficient nonviral transfection techni que for human bone marrow-derived mesenchymal stem cells[END_REF][START_REF] Ferreira | Optimization of a gene electrotransfer method for mesenchymal stem cell transfection[END_REF]. In this work we achieved cell viabilities up to 90 % and cell recoveries near 50 %. Overall we were able to obtain a yield of transfection of 30 %.

We further evaluated the ability of microporated cells expressing GFP to form colonies under antibiotic selection. The plasmid used in this study had a resistance to Zeocine which wa s added to the cell culture (250 g/ml) two days after microporation. After ten days, all non-transfected cells were dead and the first colonies were expanded maintaining around 95 % of cells continuously expressing GFP. HEK293T cells can be used either for recombinant protein production or as a model to study a specific cellular feature or disease. For both applications gene delivery might be essential and in this study we presented an approach to optimize gene delivery to these cells.

Conclusions

In conclusion, this work illustrates the feasibility of microporation for the efficient delivery of plasmids into mammalian cells, using an easy to prepare buffer, whilst maintaining cell viability concurrently with high cell recoveries. Moreover, the usefulness of CCD for gene delivery optimization was demonstrated. This approach might be successfully applied to other transfection protocols for evaluating interactions between factors involved on transfection efficiency, therefore contributing to a faster monitoring and optimization.

GFP + cells (%) 1 -1 -1 -1 -1 0.4 2 1 -1 -1 1 26.1 3 -1 1 -1 -1 2.2 4 1 1 -1 -1 0.0 5 -1 -1 1 -1 0.9 6 1 -1 1 -1 0.0 7 -1 1 1 -1 10.3 8 1 1 1 -1 0.0 9 -1 -1 -1 1 0.5 10 1 -1 -1 1

Fig. 1

 1 Fig. 1 Predicted response surface plots relating the percentage of GFP + cells with two independent

Fig. 2

 2 Fig. 2 Effect of microporation buffers SB-PBS (A) and RB (B) on gene delivery to HEK293T cells

Figure 2

 2 Figure 2

Table 1

 1 Uncoded and coded levels of the independent variables.

	Independent Variables	-1	Coded levels 0	+1
	X 1	Electric pulse (V)	900	1200	1500
	X 2	Pulse width (ms)	10	20		30
	X 3	Pulse number	1	2		3
	X 4	DNA amount ( g)	1	2		3

Table 2

 2 Experimental results of the dependent variable (GFP + cells (%)) after microporation of HEK293T using a CCD design setup.

	Run	Independent variables	Dependent variable
	X 1	*	X 2 *	X 3 *	X 4 *

  Uncoded levels of independent variables (X 1 =electric pulse; X 2 =pulse width; X 3 =pulse number; X 4 =DNA amount) are presented in Table1.

						35.2
	11	-1	1	-1	1	4.2
	12	1	1	-1	1	0.0
	13	-1	-1	1	1	1.9
	14	1	-1	1	1	17.4
	15	-1	1	1	1	8.3
	16	1	1	1	1	0.0
	17	-1	0	0	0	9.2
	18	1	0	0	0	10.1
	19	0	-1	0	0	25.0
	20	0	1	0	0	27.6
	21	0	0	-1	0	32.6
	22	0	0	1	0	29.3
	23	0	0	0	-1	34.7
	24	0	0	0	1	55.1
	25	0	0	0	0	38.6
	26	0	0	0	0	39.2
	27	0	0	0	0	35.0
	28	0	0	0	0	37.0

*

Table 3

 3 ANOVA for the reduced quadratic model (Eq. 3) that predicts the response variable GFP + cells, as a function of the independent variables after microporation.

	Source	SS	df MS F-value F-critical p-value
	Regression 7167 13 551	23	3	<0.001
	Residual	335 14 24
	Lack of fit
	(LF) 323 11 29	8	9	0.06
	Pure Error
	(PE) 12	3	4
	Total	7501 27

Acknowledgments

We would like to thank Martins, S.M.A. for helpful discussions on the CCD analysis. This work was funded by Fundação para a Ciência e a Tecnologia through Post-doctoral grants awarded to C. Madeira (SFRH / BPD / 41755 / 2007) and to S. Ribeiro (SFRH / BPD / 41824 / 2007).