A Cangiani
email: andrea.cangiani@unimib.it

R Natalini
email: roberto.natalini@cnr.it.

A spatial model of cellular molecular trafficking including active transport along microtubules

Keywords: 2000 Mathematics Subject Classification. Primary: 35K57; Secondary: 35K60, 65M60, 92C37, 92C45 Signal transduction, intracellular transport, protein motors, microtubules, Discontinuous Galerkin Methods

We consider models of Ran-driven nuclear transport of molecules such as proteins in living cells. The mathematical model presented is the first to take into account for the active transport of molecules along the cytoplasmic microtubules. All parameters entering the models are thoroughly discussed. The model is tested by numerical simulations based on Discontinuous Galerkin finite element methods. The numerical experiments are compared to the behavior observed experimentally.

Introduction

All cells receive and respond to signals from their surrounding. Any external stimulus acting on the cell plasma membrane activates several internal second messenger reactions that regulate virtually all aspects of cell behavior, including metabolism, movement, proliferation, and differentiation. As eukaryotic cells are highly compartmentalized systems in which biochemical reactions occur in physically distinct regions, the signal has to be transduced to the correct compartment were the cellular response to the external environment is initiated. All together, this process is named cellular signal transduction.

The nucleus of eukaryotic cells and specifically the genomic DNA, is the target of many intracellular transduction pathways. Indeed the response of the cell to the impinging signal is obtained through the expression of specific genes. In fact, as protein synthesis is carried out in the ribosome, the cellular response depends as well on the export of RNAs out of the nucleus.

Here, we concentrate on molecular trafficking across the nuclear envelope and nucleocytoplasmic transport. The size of the molecules involved is small enough to permit efficient diffusion in the cytoplasm and nucleus. The translocation across nuclear envelope (NE) may proceed through the nuclear pore complexes (NPCs) following two different mechanisms: passive diffusion or facilitated translocation. Only those molecules of mass less than 40 kDa [START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF] can freely diffuse through the NPCs. To permit the translocation of larger molecules, a system for active transport across the NPCs has evolved. The cargo protein equipped with a nuclear localization signal (NLS) binds to a nucleocytoplasmic transport receptor (NTR) karyopherin known as 'importin', which mediates the transport through the NE. The energy needed by facilitated translocation is provided by the Ran complex. The Ran protein is a small GTPase (see [START_REF] Quimby | The small gtpase ran: interpreting the signs[END_REF] for a review) cycling between two states: bound to guanosine triphosphate (RanGTP, active state) and to guanosine diphosphate (RanGDP, inactive state). The irreversible hydrolysis of RanGTP into RanGDP catalyzed by RanGAP in the cytoplasm maintains the Ran-cycle out of equilibrium permitting the accumulation of RanGTP in the nucleus, see [START_REF] Kopito | Reversibility in nucleocytoplasmic transport[END_REF][START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF]. This is then used to break the importin-cargo complex and thus permit the final release of cargo into the nucleus.

The mechanism of facilitated translocation permits the selective and regulated translocation of relatively large molecules. Indeed, it is understood that the combined action of importins, Ran complex and cargoes finely regulates the action of transcription factors within the nucleus.

Many mathematical models have been proposed, which qualitatively reproduce the dynamics of intracellular trafficking and signal transduction, often formulated in terms of ordinary differential equations describing the time course of the molecular concentrations (compartmental models) [START_REF] Macara | Transport into and out of the nucleus[END_REF][START_REF] Smith | Systems analysis of ran transport[END_REF][START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF][START_REF] Becskei | Quantitative models of nuclear transport[END_REF][START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF][START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF][START_REF] Riddick | The adapter importin-α provides flexible control of nuclear import at the expense of efficiency[END_REF]. These are based on schematic descriptions of signal transduction pathways which do not take into account the spatial localization of the reactions [START_REF] Harding | Subcellular localization determines map kinase signal output[END_REF]. See [START_REF] Eungdamrong | Modeling cell signaling networks[END_REF] for a review on signaling networks modeling.

Spatial simulations of cellular signal transduction pathways has been considered by some authors, see Smith et al. in [START_REF] Smith | Systems analysis of ran transport[END_REF] and the recent review by Kholodenko [START_REF] Kholodenko | Spatially distributed cell signalling[END_REF]. In [START_REF] Smith | Systems analysis of ran transport[END_REF], spatial and compartmental simulations are compared, giving similar answer. In our opinion, this has to be the case if we do not introduce any spatial details into the spatial model. In fact, if the parameters are obtained by fitting experimental data with compartmental simulations, spatial models may as well be less realistic. Let us remark here, that in the following the word spatial will be used essentially to refer to spatial gradients inside of each compartment. This is in contrast with the use of this term to refer to the gradient between the compartements.

A first attempt to build up a spatial theory of intracellular signaling was made in [START_REF] Serafini | Mathematical models for intracellular transport phenomena[END_REF], where some reaction diffusion models were proposed with a preliminary analytical discussion about local and global well-posedeness of the equations in domains with permeable membranes. In the present paper, we introduce a spatial integrated model for Ran-driven nuclear import of molecules incorporating diffusion and membrane transport for a large-scale model of living cell, by discussing the crucial problem of parameters and pathways localization. For the first time, we also take into account the active transport along the microtubules of the importin-cargo complex. The microtubules, together with the cytoplasmic filaments, constitute the cytoplasm's dynamic structure that maintains cell shape (cytoskeleton). They facilitate nuclear import of some proteins and viruses by providing a preferential way directed towards the nuclear envelope (see [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF] and references therein). This is a typically spatial phenomena that cannot be taken into account in compartmental models. Even if the actual mechanism and the scale are quite different, this investigation can be related to some models where a comparison is made of the relation between transport and diffusione, see for instance [START_REF] Loverdo | Enhanced reaction kinetics in biological cells[END_REF].

We observe also that, to simplify the model, here the NPC itself is considered to act in average as a porous membrane, with an effective permeability for each molecular species, and it is mathematically modeled by using Kedem-Katchalsky conditions and so it does not constitute an intermediate chamber in the translocation.

The system of non-linear equations arising from spatial modeling shall be solved using a new numerical technique based on Discontinuous Galerkin schemes. The details on the derivation and numerical properties of the schemes are given in [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF]. Here we shall compare our approach to previous models and perform some simulations using experimental data in a realistic framework.

Reaction-diffusion model

The present model originates from the ODE model of Ran-driven nucleocytoplasmic import successively developed by Gorlich et al. [START_REF] Ribbeck | Kinetic analysis of translocation through nuclear pore complexes[END_REF][START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF], Smith et al. [START_REF] Smith | Systems analysis of ran transport[END_REF], Riddick and Macara [START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF][START_REF] Riddick | The adapter importin-α provides flexible control of nuclear import at the expense of efficiency[END_REF], and Kopito and Elbaum [START_REF] Kopito | Reversibility in nucleocytoplasmic transport[END_REF].

Following [START_REF] Kopito | Reversibility in nucleocytoplasmic transport[END_REF], we keep the reaction network to the essentials. We simplify the set of equations and explicitly introduce the spatial component within the variables which shall represent the molar concentrations of the molecular species.

The model is a system of six coupled semilinear parabolic PDEs set on two compartments: cytoplasm and nucleus. Two equations will be added in Section 3 to take into account the translocation along microtubules, producing a system of eight coupled semilinear parabolic/hyperbolic equations.

The model includes, for each species, its initial concentration, molecular reactions, Fickian diffusion, membranes translocation conditions, and, later on, active transport along the microtubules.

In particular, facilitated translocation through the NE is only permitted to transport receptor complexes and transport along the microtubules is only permitted to complexes associated to a motor protein.

We shall make use of the following notation. We let Ω represent the cell's domain, ∂Ω its boundary (the plasma membrane), and Γ nc the interface between cytoplasm and nucleus (the nuclear envelope). Further, we let Ω c and Ω n represent the cytoplasmic and nuclear compartment, respectively. Clearly Ω n ⊆ Ω, Γ nc = ∂Ω n , and ∂Ω n ∩ Ω = ∅.

By the symbol n we will always denote the normal unit vector pointing outside of a given compartment. Finally, the jump on Γ nc of a given species concentration u will be denoted by

[[u]] := (u (c) -u (n) )| Γnc ,
where u (c) = u| Ωc and u (n) = u| Ωn .

2.1. Reaction network and mass action law. The reaction terms are written in terms of the Law of Mass Action [START_REF] Vries | A Course in Mathematical Biology: Quantitative Modeling with Mathematical & Computational Methods[END_REF]. We assume that each biochemical reaction pathway can be decomposed into unidirectional elementary reaction and that for each elementary reaction the rate of change of the reactants is proportional to the product of their concentrations. The Law of Mass Action is in fact a mathematical model expressing the fact that the reaction depends on the number of molecular collisions and the probability that a collision has enough kinetic energy to initiate the reaction. The constant of proportionality is thus named kinetic constant.

Experimental values of the kinetic constants are usually available from in vitro experimentation from purified components, and thus they do not account for environmental interactions, competition and localizations. In fact, one of the goals of spatial modeling should be to obtain the correct 'scaling relationship' of the parameters by fitting to known functional effects [START_REF] Eungdamrong | Modeling cell signaling networks[END_REF].

Another limitation to quantitative mathematical modeling regards the rates of enzymatic irreversible reactions: often the literature reports the Michaelis-Menten (MM) kinetic parameters of the catalyzed reaction instead of the kinetic constant of each reaction composing it (see [START_REF] Michaelis | Die kinetik der invertinwirkung[END_REF][START_REF] Briggs | A note on the kinetics of enzyme action[END_REF][START_REF] Segel | The quasi-steady-state assumption: A case study in perturbation[END_REF], or the review in Chapter 6 of [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF]). Again, we have to stick to the data we are given, bearing in mind that an extra approximation has been introduced in the model.

The network of reactions involved in Ran-driven nuclear import pathway are schematized in Figure 1 and described in Table 1 using the species names of Table 2.

We shall assume that facilitated translocation through the nucleocytoplasmic membrane is only permitted to molecules associated to a transport receptor (see below the description on membrane shuttling).

Receptors (adapters)-mediated import necessitates RanGTP to disassociate the cargo from the receptor once the complex has entered the nucleus. Thus, the transport mechanism relies on the presence of large concentrations of RanGTP in the nucleus in comparison to the cytoplasm (Ran gradient). The Ran complex is responsible to maintain the RanGTP/RanGDP gradient, thus permitting cargo accumulation in the nucleus.

Within the cell, the small GTPase Ran can bind to guanosine nucleotide phosphates GDP and GTP, forming the RanGDP and RanGTP complexes, respectively. These reaction, described in the model by MM-kinetics, are Rt Table 1. Kinetic constants for the reactions involved in the Ran-driven transport process, with respect to the species names in Table 2. Second order constants are measured in μM -1 s -1 while first order constants are measured in s -1 . RCC1 catalyzed reaction is represented with a MM scheme in terms of the constant concentration c1, as in Smith e al, 2002. For a more refined scheme (a multistep scheme) see [START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF]. The MM scheme gives a good approximation as the intermediate steps (complex formation) are more rapid than the exchange reaction. The enzymes concentrations are assumed to remain constant. In particular, following [START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF], we set R g = 0.5μM and C 1 = 0.7μM.

Rg → Rd Ω c m1(Rt) := K t c Rg Rt K t m + Rt K t c K t m 10.6 s -1 0.7 μM [23] Rd C 1 → Rt Ω n m2(Rd) := K d c C1 Rd K d m + Rd
catalyzed by two specific enzymes: RanGAP which is located in the cytoplasm and RanGEF (RCC1) located within the nucleus. Due to this cycle of reactions, a large concentration jump of RanGTP across the NE is readily established, with a high RanGTP concentration in the nucleus. The first cycle of cargo import begins with the formation of the cargo receptor complex. The receptor involved in nuclear import is made of two subunits: the adaptor Importin-α which binds the cargo Nuclear Localization Signal (NLS) and the receptor Importin-β which permits the translocation through the NPC. Here, we simplify matters by considering the binding of the cargo with the whole importin complex (see [START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF][START_REF] Riddick | The adapter importin-α provides flexible control of nuclear import at the expense of efficiency[END_REF] for a detailed analysis of the various import pathways).

The second (parallel) cycle is the nuclear import of cytoplasmic RanGDP by the effector NTF2. Here we assume that NTF2 is always available and treat RanGDP as all bound to NTF2. Nuclear RanGDP, imported by the transport receptor NTF2, interacts with the nuclear RanGEF (RCC1) which catalyzes the exchange between nucleotide GDP and GTP forming the nuclear RanGTP. The complete biochemical scheme of this exchange is composed of a sequential set of enzymatic steps [START_REF] Klebe | The kinetic mechanism of rannucleotide exchange catalyzed by rcc1[END_REF][START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF][START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF] which, following [START_REF] Smith | Systems analysis of ran transport[END_REF], we schematize as a single reaction.

The nuclear importin-cargo complex can bind the RanGTP to form two complexes, RanGTP•Impβ and Impα•cargo. This latter complex can dissociate, create free Impα and cargo or interact with RanGTP to form the RanGTP•Impα complex and free cargo. We simplify this pathway with the single reaction

T • C + R → C + T • R,
written in terms of the generic transport receptor T , cargo C, and with R representing the concentration of Ran (bound to GTP in this case). The free cargo, which usually is an activated transcription factor, can now bind DNA and activate the gene expression program. In our model, these reactions are neglected, and thus the cargo is allowed to accumulate within the nucleus, while the receptor-RanGTP complex exits the nucleus and, eventually, dissociates.

Diffusion. The molecular species diffusion is expressed in terms of

Fick's law. All the molecular species diffuse with a specific diffusion coefficient. This is obtained for low Reynolds number from the Einstein relation

d = KT 6πηR s ,
in terms of the Boltzmann constant K, the absolute temperature T , the Stokes radius R s , and the viscosity of the medium η (see e.g. [START_REF] Ribbeck | Kinetic analysis of translocation through nuclear pore complexes[END_REF]).

The viscosity of the cytoplasm has been measured in the nineteen-nineties for a wide range of molecular weights and cellular environments [START_REF] Seksek | Translational diffusion of macromolecule-sized solutes in cytoplasm and nucleus[END_REF]. It was initially thought that the viscosity had to depend on the molecular size to account for the sieving effect due to the cytoskeletal network and other macromolecular structures. For this reason, solutes diffusion was described by the translational diffusion coefficient [START_REF] Kao | Determinants of the translational mobility of a small solute in cell cytoplasm[END_REF], depending on the viscosity of the fluid-phase cytoplasm and on the collisions with intracellular components. The former, representing the viscosity sensed by a small solute in the absence of interactions with macromolecules and organelles, is 1.2-1.4 times that of water [START_REF] Fushimi | Low viscosity in the aqueous domain of cell cytoplasm measured by picosecond polarization microfluorimetry[END_REF].

The net viscosity of the cytoplasm of Swiss 3T3 fibroblasts with respect to small solutes such as metabolites and nucleic acids was measured by Kao et al. [START_REF] Kao | Determinants of the translational mobility of a small solute in cell cytoplasm[END_REF] using fluorescence recovery after spot-photobleaching (FRAP) [START_REF] Axelrod | Mobility measurement by analysis of fluorescence photobleaching recovery kinetics[END_REF]. They found that mobility was around four to fivefold slower than in water. Contrary to expectations, later studies extended the validity of this result to the mobility of macro-molecules of up to few hundreds of kDs of mass [START_REF] Seksek | Translational diffusion of macromolecule-sized solutes in cytoplasm and nucleus[END_REF]. Moreover, the same applies to the nucleus [START_REF] Seksek | Translational diffusion of macromolecule-sized solutes in cytoplasm and nucleus[END_REF] and the mitochondrial matrix [START_REF] Partikian | Rapid diffusion of green fluorescent protein in the mitochondrial matrix[END_REF], thus viscosity fivefold higher than water can be assumed throughout.

A correction may be employed in the proximity of the plasma membrane as the translational diffusion of the fluorescent probe BCECF already used by [START_REF] Kao | Determinants of the translational mobility of a small solute in cell cytoplasm[END_REF] is found to be twofold lower near the cell membrane due to high density of proteins [START_REF] Swaminathan | Cytoplasmic viscosity near the cell plasma membrane: translational diffusion of a small fluorescent solute measured by total internal reflection-fluorescence photobleaching recovery[END_REF].

The diffusion coefficients used in the model are reported in Table 2. 2. Diffusion coefficients for the molecular species involved in the Ran-driven transport system. Values are the same for the cytoplasm and the nucleus. These data have been calculated assuming a cytoplasm viscosity 5 times higher than water [START_REF] Seksek | Translational diffusion of macromolecule-sized solutes in cytoplasm and nucleus[END_REF] at a temperature of 20 • Celsius. Cargo weight taken as the ERK mass, 42 kDa (for ERK1) or 44 kDa (ERK2). Data from the literature are not always consistent: for instance, Smith et al [START_REF] Smith | Systems analysis of ran transport[END_REF] estimates the diffusivity of Ran in 30μm 2 /s.

Diffusion coefficients for

Facilitated translocation through the nuclear envelope.

In this paper, we shall focus on nuclear import by assuming that the cargo is already in the cell and that there is no exchange of substances between the cell and the surrounding environment through the plasma membrane. Mathematically, no-flux conditions are specified for all species at the plasma membrane. Thus, given the generic species concentration u, we impose the following condition at the plasma membrane:

(1)

d u ∂u ∂n | ∂Ω = 0,
where, as always in these pages, n represents the unit normal vector pointing outside the cell. In fact, the mechanism of transport through the plasma membrane and the NE are quite different. Most often, external signals are transduced from external membrane receptors into internal second messenger activation on the inside of the plasma membrane [START_REF] Alberts | Molecular Biology of the Cell[END_REF]. Thus, there is no passage of matter through the membrane. On the contrary, here we are interested in signals transduced across the nuclear membrane by a mass transfer mechanism. All shuttling across the nuclear envelope takes place through the nuclear pore complexes (NPCs), multi-component protein structures spanning the nuclear envelope. Two mechanisms of translocation are permitted: passive diffusion or facilitated translocation. Due to the limited diameter of the pore lumen (about 10 nm), ions and small metabolites can freely diffuse through the NPCs only if their mass is less than ∼ 40 kDa. Molecules of size greater than that of the NPCs can still shuttle across the NE by facilitated translocation. The actual mechanism of facilitated translocation is still under research (for a review of various models, see [START_REF] Fahrenkrog | The nuclear pore complex: nucleocytoplasmic transport and beyond[END_REF]). It is, though, well established that the molecule must posses a specific aminoacid domain which binds to transport receptors (importins and exportins) that operate a conformational change in the NPC whose 'functional size' is in the order of 25 nm [START_REF] Mattaj | Nucleocytoplasmic transport: The soluble phase[END_REF][START_REF] Becskei | Quantitative models of nuclear transport[END_REF]. Thus, we assume that facilitated translocation is allowed only to receptor complexes [START_REF] Kopito | Reversibility in nucleocytoplasmic transport[END_REF].

The study presented in [START_REF] Kopito | Reversibility in nucleocytoplasmic transport[END_REF][START_REF] Kopito | Nucleocytoplasmic transport: a thermodynamic mechanism[END_REF] proves that the accumulation of cargo in the nucleus is not due, as previously thought, to the ability of the 'importin' receptor to cross the NE being unidirectional (from the cytoplasm to the nucleus). The translocation is actually bidirectional, with the accumulation of cargo in the nucleus being rather due to the asymmetric concentration of RanGTP, as explained above.

In the absence of a detailed and realistic model of the reactions happening within the pore (interactions of proteins with nucleoporins, meshwork of filaments within the lumen of the pore, etc), a more convenient approach is to use a 'coarse-grain' formulation, in terms of permeabilities times concentration differences between the two compartments (see [START_REF] Smith | Systems analysis of ran transport[END_REF]). The flux across the nucleocytoplasmic boundary is modeled as the product of a proportionality constant (the permeability P ) times the concentration difference across the nucleocytoplasmic boundary, as a direct consequence of the fact that the flux does not require additional energy input [START_REF] Smith | Systems analysis of ran transport[END_REF]. Thus, we fix the following Kedem-Katchalsky transmission conditions [START_REF] Kedem | Thermodynamic analysis of the permeability of biological membrane to non-electrolytes[END_REF] (see also [START_REF] Quarteroni | Mathematical and numerical modeling of solute dynamics in blood flow and arterial walls[END_REF][START_REF] Calabrò | Global existence of strong solutions of parabolic problems on partitioned domains with nonlinear conditions at the interface. applications to mass transfer processes through semi-permeable membranes[END_REF][START_REF] Serafini | Mathematical models for intracellular transport phenomena[END_REF]] for a mathematical and numerical treatment of this conditions) at the nuclear envelope:

(2)

d u ∂u (n) ∂n | Γnc = p u [[u]],
where here the normal vector n is pointing from the nucleus to the cytoplasm. The nuclear membrane permeability is assumed constant as experimental data indicate binding to the pore complex is far from saturation [START_REF] Ribbeck | Kinetic analysis of translocation through nuclear pore complexes[END_REF][START_REF] Smith | Systems analysis of ran transport[END_REF]. The permeability values can be calculated from experimental values on the capacity, area, and number of the NPCs present on the nuclear envelope. Ribbeck and Gorlich [START_REF] Ribbeck | Kinetic analysis of translocation through nuclear pore complexes[END_REF] estimated > 100 translocation events/NPC/second and a density of NPCs of 5.1 ± 0.2 NPCs/μm 2 in the NE of HeLa cells. Estimates of translocation rates of most molecular species are available from the literature and have already been used in a number of studies [START_REF] Macara | Transport into and out of the nucleus[END_REF][START_REF] Smith | Systems analysis of ran transport[END_REF][START_REF] Görlich | Characterization of ran-driven cargo transport and the rangtpase system by kinetic measurements and computer simulation[END_REF][START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF][START_REF] Riddick | The adapter importin-α provides flexible control of nuclear import at the expense of efficiency[END_REF]. The values used in the model are taken from [START_REF] Smith | Systems analysis of ran transport[END_REF] and are shown in Table 3. These were fitted by comparison of experimental data with compartmental simulations: more investigation on the correct values for spatial modeling is needed, for instance by following the approach of the homogenized 'effective' permeabilities proposed in [START_REF] Berezhkovskii | Boundary homogenization for trapping by patchy surfaces[END_REF]. Table 3. Permeability values for molecular species involved in membrane translocation. The permeability is the kinetic effect of the molecular species with the nucleoporin complex (NPC).

The set of transmission conditions at the NE is closed by imposing continuity of the flux:

(3)

d u ∂u (c) ∂n = d u ∂u (n) ∂n .
Here, we use the same normal vector for both sides to emphasize the identity of the fluxes.

For molecular species which cannot pass through the nuclear membrane we just impose homogeneous Neumann boundary conditions on the relevant compartment. For instance, if the species u is confined in the cytoplasm, then (4)

d u ∂u (c) ∂n = 0.
2.4. Transport model for diffusive molecules. Most signaling cargoes resort on diffusion only to reach their target cellular compartment [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF]. Thus, before undertaking into the inclusion of active transport, it is relevant to display the basic diffusion model of Ran-driven nuclear import. This is based on the building elements discussed up to now, namely diffusion, complex-forming reactions, and active NE translocation. A similar spatial diffusion model of Ran-driven nuclear import has already been presented in [START_REF] Smith | Systems analysis of ran transport[END_REF].

By collecting all the model's terms described above, we obtain the following system of coupled semilinear parabolic equations.

In the cytoplasmic compartment Ω c , the species concentrations obey:

(5)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂R t ∂t = d r ΔR t -m 1 (R t ) -r 1 (R t , T ) + r -1 (T r ), ∂R d ∂t = d r ΔR d + m 1 (R t ), ∂T r ∂t = d tr ΔT r + r 1 (R t , T ) -r -1 (T r ), ∂C ∂t = d c ΔC -r 2 (C, T ), ∂T ∂t = d t ΔT -r 1 (R t , T ) + r -1 (T r ) -r 2 (C, T ), ∂T c ∂t = d tc ΔT c + r 2 (C, T ).
We refer to Tables 1 and2, for the notation about concentrations and for the values of the various constants. As we assume that no matter is entering or exiting the cell through the plasma membrane, the homogeneous Neumann boundary condition ( 1) is imposed along ∂Ω to all species.

In the nuclear compartment Ω n , we have the following system of reactiondiffusion equations:

(6) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂R t ∂t = d r ΔR t + m 2 (R d ) -r 1 (R t , T ) + r -1 (T r ) -r 3 (R t , T c ), ∂R d ∂t = d r ΔR d -m 2 (R d ), ∂T r ∂t = d tr ΔT r + r 1 (R t , T ) -r -1 (T r ) + r 3 (R t , T c ), ∂C ∂t = d c ΔC + r 3 (R t , T c ), ∂T ∂t = d t ΔT -r 1 (R t , T ) + r -1 (T r ), ∂T c ∂t = d tc ΔT c -r 3 (R t , T c ),
The two systems of equations above are coupled through the appropriate transmission conditions on Γ nc :

(7) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d r ∂R (c),(n) t ∂n = 0, -d r ∂R (c) d ∂n = d r ∂R (n) d ∂n = p d [[R d ]], -d tr ∂T (c) r ∂n = d tr ∂T (n) r ∂n = p tr [[T r ]], d c ∂C (c),(n) ∂n = 0, -d t ∂T (c) ∂n = d t ∂T (n) ∂n = p t [[T ]], -d tc ∂T (c) c ∂n = d tc ∂T (n) c ∂n = p tc [[T c ]].
Notice that now we are consistent with our notations and the normal vector n is relative to the two different compartements, cytoplasm for the first terms and nucleus for the second. These equations need to be provided with initial conditions which are specified below in Section 5.

Active transport along the microtubules

The microtubules (MTs) are hollow cylindrical filaments (internal and external diameter of approximately 15 and 25nm, respectively) which, together with the cytoplasmic filaments, constitute the cytoskeleton. During most of the life of the cell (interphase), the MTs are organized within the cytoplasm as an aster originating from a microtubule organizing center (MTOC) located in the proximity of the nucleus. Further, they are characterized by an orientation, with the plus ends located near the MTOC.

One of the many roles played by the MTs is that of enhancing intracellular trafficking. Indeed, the MTs are known to be used as preferential ways of motion by macro molecules (e.g. adenoviruses [START_REF] Dinh | Theory of spatial patterns of intracellular organelles[END_REF]) and intracellular organelles (e.g. endocytic vesicles [START_REF] Dinh | Theory of spatial patterns of intracellular organelles[END_REF]). The size of these entities (the adenoviruses are 90nm in diameter) limits considerably their diffusion speed in the cytoplasm; for this reason they must resort to active transport along the MTs in order to reach their target location.

It is by now well established that also some small molecules utilize motorassisted transport along the MTs [START_REF] Gundersen | Microtubules and signal transduction[END_REF][START_REF] Campbell | Role of the cytoskeleton in nuclear import[END_REF][START_REF] Salman | Nuclear localization signal peptides induce molecular delivery along microtubules[END_REF]. Examples goes from mRNA [START_REF] Tekotte | Intracellular mrna localization: motors move messages[END_REF] localizing after nuclear export, to the tumor suppressor proteins p53 [START_REF] Giannakakou | p53 is associated with cellular microtubules and is transported to the nucleus by dynein[END_REF] and Rb [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF], and the Parathyroid hormone-related protein PTHrP [START_REF] Lam | Nuclear Transport of Parathyroid Hormone (PTH)-Related Protein Is Dependent on Microtubules[END_REF]. We remark that many similarly sized molecules resort on passive diffusion only, thus, in this respect, active transport is not essential to traffiking precesses. Rather it must be seen as a way to improve its efficiency. Our aim is to introduce a complete model of transport mechanisms, accounting for both diffusion and active transport along the micrutubules, that can be used to better understand the characteristics and importance of the latter mechanism.

Active transport along the MTs is permitted by binding to a motor protein [START_REF] Nédélec | Dynamic concentration of motors in microtubule arrays[END_REF][START_REF] Smith | Models of motor-assisted transport of intracellular particles[END_REF], which possesses a mechanism for moving along the MT at a speed of about 0.5 to 1 μm s -1 [START_REF] Smith | Models of motor-assisted transport of intracellular particles[END_REF][START_REF] Nédélec | Dynamic concentration of motors in microtubule arrays[END_REF]. Two families of motor proteins associate to the MTs: dynein, which permits transport from the plus end to the minus end, and kinesin, which transports in the opposite direction. In the case of viruses and intracellular organelles, changes of direction of motion and pauses are often observed, suggesting that molecules may proceed by repeated detaching and changing type of motor [START_REF] Smith | Models of motor-assisted transport of intracellular particles[END_REF]. On the contrary, proteins such as the tumor suppressor protein p53 are known to bind to the MTs in association with the motor protein dynein only [START_REF] Giannakakou | Enhanced microtubule-dependent trafficking and p53 nuclear accumulation by suppression of microtubule dynamics[END_REF]. Thus we can assume that, once associated to a motor protein and until detachment, a cargo protein will move towards the MTOC at a fixed speed of about 1 μm s -1 . The actual mechanism of association to an MT is still under scrutiny, with tentative mathematical models to fully describe MTs effective directional transport being proposed, for instance, in [START_REF] Smith | Models of motor-assisted transport of intracellular particles[END_REF][START_REF] Nédélec | Dynamic concentration of motors in microtubule arrays[END_REF][START_REF] Dinh | A model for intracellular trafficking of adenoviral vectors[END_REF][START_REF] Salman | Nuclear localization signal peptides induce molecular delivery along microtubules[END_REF][START_REF] Kahana | Active transport on disordered microtubule networks: The generalized random velocity model[END_REF]. Here we adopt the following continuous model which will be validated by comparison with the experimental results of Roth et al. [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF] and Lam et al. [START_REF] Lam | Nuclear Transport of Parathyroid Hormone (PTH)-Related Protein Is Dependent on Microtubules[END_REF].

We introduce two new species concentrations: the dynein D and the dynein, cargo, and importin complex D c . The dynein is free to diffuse in the cytoplasm (with speed comparable to that of Ran [START_REF] Nédélec | Dynamic concentration of motors in microtubule arrays[END_REF]), and thus its dynamics will be described by a reaction-diffusion equation similar to those previously considered. A cargo and importin complex can associate to a dynein molecule to form the complex D c and in this way acquiring the ability to be transported along an MT. We actually identify attachment/detachment to the MT with the reaction forming/breaking the complex D c , which are given by the kinetic reactions described in Table 4. In other words, we assume that the complex D c is only found bound to the MTs, and we treat it by an hyperbolic equation of transport coupled to the other, parabolic, equations through the reactions of Table 4. Table 4. The reactions of attachment to and detachment from the microtubules are identified with the forming and breaking of the complex D c .

Reactions and kinetic constants for the dynein complex

Let b : Ω c → R 3 be a velocity field associated with the active transport along the MTs. We introduce the following differential equations modeling

Ω c Ω c Ω c Ω n
⎧ ⎪ ⎨ ⎪ ⎩ ∂D ∂t = d d ΔD -r a (T c , D) + r d (D c ), ∂D c ∂t = -∇(b D c ) + r a (T c , D) -r d (D c ). (8)
Accordingly, we modify the equation for the concentration T c of the importincargo complex in the cytoplasm, by adding the MTs related reactions:

(9) ∂T c ∂t = d tc ΔT c + r 2 (C, T ) -r a (T c , D) + r d (D c ),
It remains to specify the nature of the advective field b and, accordingly, the boundary conditions for the new species. For simplicity, we assume that b is a radial field pointing to a point inside the nucleus, cf. Figure 2. To reproduce the fact that the MTs are not present up to the nuclear envelope, we let Ω c ⊂ Ω c denote a portion internal to the cytoplasm where the field b is of constant modulus equal to 1, and we assume that on the remaining layer Ω c = Ω c \Ω c the field b drops to zero with continuity, so that b| Γnc = 0. We make the layer region Ω c coincide with the end points of the MTs by imposing detachment. This can be done assuming that the attachment rate r a drops to zero in Ω c , while the detachment rate T d goes to infinity towards Γ nc (how these condition is handled in practice will be described in the section on numerical simulation).

The new species D and D c are present only on Ω c , thus we need to provide boundary conditions for both ∂Ω and Γ nc . We start by noticing that the boundary Γ nc can be seen as an outflow boundary for b (actually, b| Γnc = 0 and our reactions rates imply that D c = 0 along Γ nc ). Thus we do not need to impose any boundary condition for D c along Γ nc . All other situations are to be treated with homogeneous Neumann boundary conditions, namely:

(10) ⎧ ⎨ ⎩ d d ∂D (c) ∂n | Γnc∪∂Ω = 0 ∂D (c) c
∂n | ∂Ω = 0 Finally, let us introduce the initial conditions which complete the mathematical model. The initial conditions are the concentrations of the single molecular species when cell is at rest (i.e. in the absence of external stimuli). In such a situation, the concentrations of the complexes are zero (i.e. no complex is formed before the stimulus activates the cargo), and the experimentally observed RanGTP accumulation of RanGTP in the nucleus forms only during the initial phase of the experiment. The initial concentrations used are reported in Table 5. Various experiments based on different NLS cargo injections will be considered, thus the initial concentration of the NLS cargo will be discussed in Section 5.

Initial concentrations reactant localization μM reference(s)

R t cyto 3 [49, 56] R d cyto 3 [49, 56] T cyto 4 [49] D cyto 1 n/a
Table 5. Initial concentrations of constituents of the Ran system used in the model by Riddick et al. [START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF]. All other initial concentrations are set to zero. The concentrations of the enzymes RCC1 and RanGAP are assumed to be constant, as described in Section 2.1.

Numerical approximation of the mathematical model

In order to ease notation in the presentation of the numerical method, we rewrite the mathematical model given by equations ( 5), ( 6), ( 7), ( 8), [START_REF] Briggs | A note on the kinetics of enzyme action[END_REF], and (10) in vector form.

The pde problem in vector form. Let us denote the number of unknown concentrations by n, and collect all concentrations variables in the vector function

u := (u 1 , . . . , u n ) T : Ω c ∪ Ω n → R n .
We shall look for solutions belonging to the space

H 1 := [H 1 (Ω c ∪ Ω n )] n at any time t ∈ [0, T ],
with T representing some final time.

We introduce the diagonal matrix U = diag(u 1 , . . . , u n ) and the gradient ∇u : Ω c ∪ Ω n → R n×d given by ∇u := (∇u 1 , . . . , ∇u n ) T , with

∇u i (x) ∈ R d , i = 1, . . . , n, x ∈ Ω c ∪ Ω n . Further, for a tensor Q : Ω c ∪ Ω n → R n×d , with rows Q i , i = 1, . . . , n, we define its divergence ∇ • Q : Ω c ∪ Ω n → R n by ∇ • Q := (∇ • Q 1 , . . . , ∇ • Q n ) T .
Our model couples diffusion equations such as [START_REF] Briggs | A note on the kinetics of enzyme action[END_REF] with the transport equation appearing in [START_REF] Berezhkovskii | Boundary homogenization for trapping by patchy surfaces[END_REF], which is only defined in the cytoplasm Ω c . In order to write all equations at once, we collect the advection and diffusion coefficients in, respectively, the tensor B ∈ [C 1 (0, T ; Ω)] n×d , whose rows are denoted by B i , i = 1, . . . , n, and the diagonal tensor

A ∈ [C(0, T ; Ω c ∪ Ω n )] n×n with A = diag(d 1 , d 2 , . . . , d n )
, where d i ≥ 0, i = 1, . . . , n are the diffusivity of the various species. In accordance to our model, we assume that B| Ωn ≡ 0 and Bn ≤ 0. Further, we define the diagonal tensor P := diag(p 1 , . . . , p n ) where p i ≥ 0, i = 1, . . . , n are the species permeabilities. It is assumed that, if a given term does not apply to a given species, the relevant coefficient is set to zero. For instance, B i will be identically zero if the species u i is not transported by the MTs.

Finally, we collect all reaction terms in the vector field r(u), which depends on both the species and the compartment, and let u 0 ∈ [L 2 (Ω)] n represent the initial conditions.

We look for solutions u ∈ L 2 (0, T ; H 1 ) of the following initial and boundary value problem:

(11) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u t = ∇ • (A∇u -U B) + r(u) in [0, T ] × (Ω c ∪ Ω n ), u(0, x) = u 0 (x) o n{0} × Ω, A∇u -Bu n = 0 on ∂Ω,
A∇u n| Ωc = P (u| Ωn -u| Ωc ) on Γ nc , A∇u n| Ωn = P (u| Ωc -u| Ωn ) on Γ nc .

More general boundary conditions, including mixed-type conditions, and transmission conditions, including nonlinear conditions, are discussed in [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF].

The second equation in [START_REF] Berezhkovskii | Boundary homogenization for trapping by patchy surfaces[END_REF] describing the dynamics of the dynein, cargo, and importin complex is purely hyperbolic, thus corresponding to one of the equations in [START_REF] Campbell | Role of the cytoskeleton in nuclear import[END_REF] having zero diffusion coefficient. As mentioned above, such equation is solved on the whole cytoplasmic domain. To model detachment from the MTs in the proximity of the nuclear envelope, i.e. in the region Ω c , we correct the coefficients T d , r a , and b. Numerically, we cannot let T d approach infinity, as this would make for an arbitrarily stiff problem. We rather tune the restriction of T d , r a , and b to Ω c in accordance with the mesh size in order to ensure that the method is stable and D c | Γnc is as close as possible to zero. Further, to ensure conservation of mass, we correct the transmission condition for T c and D as follows:

-d tc ∂T (c) c ∂n | Γnc = d tc ∂T (n) c ∂n | Γnc = p tc [[T c ]] -p tc D c c | Γnc , -d d ∂D (c) c ∂n | Γnc = -p tc D c c | Γnc , with the extra term -p tc D c
c | Γnc to be included in the forcing term r(u).

The interior penalty discontinuous Galerkin method.

In recent years, there has been an increasing interest in discontinuous Galerkin finite element methods (DGFEM) [START_REF]Discontinuous Galerkin Methods[END_REF]. These were first introduced in 1973 by Reed and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] for the solution of hyperbolic equations, and later generalised to the solution of elliptic and parabolic pdes under the name of interior penalty (IP) methods [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF].

The renewed interest in these methods is due to their very good stability properties when used to approximate solutions to convection-dominated convection-diffusion problems, as well as due to the great flexibility in grid design they offer. Moreover, they naturally embed good local conservation properties of the state variable, which can be advantageous for timedependent problems.

DGFEMs are suited to the solution of heterogeneous problems compraising both parabolic and hyperbolic equations. Moreover, the weak imposition of boundary conditions typical of discontinuous methods permits a very natural imposition of the transmission conditions through the choice of numerical fluxes.

The method used here, which is detailed in [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF], is a generalization of the symmetric version of the interior penalty discontinuous Galerkin (IPDG) method [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF]. This was analysed in [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF] for the solution of second-order partial differential equations with nonnegative characteristic form and later extended to the solution of semilinear time-dependent problems in [START_REF] Lasis | hp$-version discontinuous galerkin finite element method for semilinear parabolic problems[END_REF]. In view of the present application, the IPDG method is extended in [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF] to the solution of systems of equations of changing type and with transmission conditions.

The finite element method is introduced by defining a shape-regular subdivision T of Ω c ∪ Ω n into disjoint open elements κ ∈ T. We assume that the elements of T 'belong' to either subdomain, so that subdivisions of Ω c and Ω n are authomatically introduced. Abusing notation, we use the simbol Γ nc to denote the union of elemental faces belonging to both subdivisions. Further we decompose the subdivision skeleton Γ := ∪ κ∈T ∂κ into three disjoint subsets

Γ = ∂Ω ∪ Γ int ∪ Γ nc ,
where Γ int := Γ\(∂Ω ∪ Γ nc ).

We assume that the subdivision T is constructed via mappings F κ , where F κ : κ → κ are smooth maps with non-singular Jacobian, and κ is the reference d-dimensional simplex or the reference d-dimensional (hyper)cube; the maps are assumed to be constructed so as to ensure that the union of the closures of the elements κ ∈ T forms a covering of the closure of Ω, i.e., Ω = ∪ κ∈T κ.

For a nonnegative integer m, we denote by P m (κ), the set of polynomials of total degree at most m if κ is the reference simplex, or the set of all tensor-product polynomials on κ of degree k in each variable, if κ is the reference hypercube. We consider the hp-discontinuous finite element space [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF] S := {v ∈ L 2 (Ω) : v| κ • F κ ∈ P mκ (κ), κ ∈ T}.

Next, we introduce some trace operators. Let κ + , κ -be two (generic) elements sharing an edge e := ∂κ + ∩ ∂κ -⊂ Γ int ∪ Γ nc . Define the outward normal unit vectors n + and n -on e corresponding to ∂κ + and ∂κ -, respectively. For functions q : Ω → R n and Q : Ω → R n×d that may be discontinuous across Γ, we define the following quantities. For q + := q| e⊂∂κ + , q -:= q| e⊂∂κ -, and

Q + := Q| e⊂∂κ + , Q -:= Q| e⊂∂κ -, we set {q} := 1 2 (q + + q -), {Q} := 1 2 (Q + + Q -),
and

[[[q]]] := q + ⊗ n + + q -⊗ n -, [[Q]] := Q + n + + Q -n -,
where ⊗ denotes the standard tensor product operator, whereby q ⊗ w = qw T . If e ∈ ∂κ ∩ Γ ∂ , these definitions are modified as follows

{q} := q + , {Q} := Q + , [[[q]]] := q + ⊗ n, [[Q]] := Q + n.

Finally, we introduce the mesh quantities

h : Ω → R, m : Ω → R, by h(x) = diam κ, m(x) = m κ , if x ∈ κ, and h(x) = {h}, if x ∈ Γ, m(x) = {m}, if x ∈ Γ, respectively.
In order to define the IPDG finite element method, we further introduce Σ := diag(σ 1 , . . . , σ n ) the diagonal matrix containing the discontinuitypenalisation parameters σ i : Γ\Γ nc , and denote

B := 1/2 diag(|B 1 •n|, . . . , |B n • n|).
The IPDG-in-space method for the system (11) reads: find

u h ∈ L 2 (0, T ; [S] n ) such that (13) (u h ) t , v h +B(u h , v h ) = f (u h ), v h , for all v h ∈ L 2 (0, T ; [S] n ),
where •, • is the L 2 scalar product and the bilinear form B(•, •) is given by ( 14)

B(u h , v h ) := κ∈T κ (A∇u h -U h B) : ∇v h + Γnc P[[[u h ]]] : [[[v h ]]] - Γ int {A∇u h -U h B} : [[[v h ]]] + {A∇v h } : [[[u h ]]] + Γ int (Σ + B)[[[u h ]]] : [[[v h ]]] + ∂Ω (U h B) : [[[v h ]]].
A fully discrete formulation is obtained from (13) by using any standard time stepping method like, for instance any Runge-Kutta time stepping. Full details on the derivation and analysis of the method can be found in [START_REF] Cangiani | Discontinuous galerkin methods for mass transfer through semi-permeable membranes[END_REF].

For our numerical experiments we employed second order DIRK time stepping with Δt = 0.1 and linear finite elements for the space discretization. This setting ensures stable solution of the given problem. The resulting IPDG finite element algorithm has been coded using the C++ library deal.II [START_REF] Bangerth | II Differential Equations Analysis Library[END_REF][START_REF] Bangerth | deal.II -a general-purpose objectoriented finite element library[END_REF].

Numerical simulations

In Section 2 we introduced a basic reaction-diffusion model of cellular molecular trafficking which we expanded in Section 3 to incorporate a model of active transport along the MTs.

The basic model can be used to simulate nuclear import of the majority of molecules which do not bind to motor proteins [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF]. The same model, when applied to those molecular species that do bind to motor proteins, can be used to simulate the inhibition of MTs binding. By comparing the in silico results obtained by the two models, we are able to quantify the net contribution of active transport.

A technique for binding inhibition in vivo or in vitro is not known and the experimentation in the laboratory of the effectiveness of MTs active transport is based on much more invasive techniques.

For instance, the experimental results presented by Roth et al. [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF] for the protein Rb and by Lamet al. [START_REF] Lam | Nuclear Transport of Parathyroid Hormone (PTH)-Related Protein Is Dependent on Microtubules[END_REF] for the protein PTHrP are based on MT-depolymerization. Two experiments are presented in [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF]. In the first, the Rb cargo protein, which is known to bind to the MTs, is activated in a cell with intact cytoskeleton, and nuclear accumulation is recorded in time using fluorescence techniques. In a second experiment, the cell is treated with the MT-depolymerizing agent nocodazole (NCZ) to create a MT-less environment, and thus obtaining a measure of nuclear accumulation in the absence of active transport. This experiment is invasive, as MTs depolymerization changes substantially the cell environment. The full consequences of depolymerization on the cell physical environment and physiology are unknown, and thus comparison with the non invasive in silico experimentation is crucial. Indeed, as descrbed above, in silico we can easily inhibit MTs binding by turning off the corresponding reaction in the model with no consequences whatsoever on the model cell environment, apart from not permitting the association to the MTs.

Here we present the in silico results of the with active transport and inhibited binding experiments for the tumor suppressor protein Rb.

The protein Rb weights 110kDa [START_REF] Lohmann | Rb1 (retinoblastoma)[END_REF]: estimating a Stokes radius of ∼81 Å gives a diffusion coefficient of 12 μm 2 /s.

In our first experiment, we assume that Rb is uniformly distributed in the cytoplasm before activation, fixing the initial value C c (x, 0) = 0.5.

The in silico concentration of the Rb cargo at fixed times after activation is shown in Figure 3.

The accumulation in time of the protein Rb in the nucleus after activation is shown in Figure 4. Enhanced accumulation rate is evident and, et least qualitatively, confirms the behavior experimentally obtained in [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF]. In fact, in the MTs inhibited binding case, we observe higher nuclear accumulation than in the corresponing experiment of [START_REF] Roth | A microtubulefacilitated nuclear import pathway for cancer regulatory proteins[END_REF]. A quantitative comparison of in silico and in vivo results would permit to better quantify the consequences of depolymerization. The results shown were obtained on the 2D-cell en- vironment depicted in Figure 2. We have verified that these results are representative of a realistic cell environment by performing the same experiment on a (computationally much heavier) 3D-cell environment simulation. The qualitative behavior is confirmed (data not shown). Figure 5 shows the concentration after 100s of some of the molecular species involved in the 3D-cell simulation.

In our second experiment, we assume that the protein Rb is activated in a peripheral zone of the cytoplasm, as shown in the first snapshots on the left of Figure 6. The Rb concentration at subsequent times is shown in the same figure: notice in particular the accumulation towards the MTs aster after 50s from activation. To better appreciate the effectiveness of the transport mechanism, the accumulation in time in the nucleus is plotted in Figure 7.

A further advantage of in silico experimentation is that we can display the concentration of all substances involved in the transport mechanism. For instance, in Figure 8 we are able to show in two separated plots the concentration of the cargo importin complex T c and that of the cargo attached to the MTs, namely, D c after few seconds from the activation of the cargo: observe that, although D c obeys a purely hyperbolic equation, the solution is smoothed due to the coupling with the diffusive T c through the reaction terms.

Finally in Figure 9 we show the evolution of the concentration of T c in a 3D-simulation of the peripheral cargo activation for the first 40 seconds. The inhomogeneous initial condition enhance the difference between the case with and without transport, by contrast with the homogeneous situation shown in Figure 5.

Conclusions

In this paper we have presented a spatial model for intracellular molecular traffiking. This model takes into account not only the facilitated translocation through the nucleus envelope, due to the action of the importin molecules, but also the work done by proteic motors along the microtubules network, at least for some families of proteins. Actually, the novel tool given by our advection-diffusion model will permit in the future to examine and verify new experimental results. The results shown confirms that the efficiency of nuclear transport is enhanced by MTs active transport, even if the improvement is somewhat obfuscated by the high cytoplasmic diffusivity values we found in the literature. These data need confirmation due to the mentioned limits of the laboratory experiments. Having said this, a more thorough analysis of the consequences of depolymerization and of the contribution of the advection phenomenon would need a direct comparison of in silico and in vivo results. The main weak point from the experimental viewpoint is in fact in the lack of reliable data. For instance, the model seems to predict that the diffusivity used is over-estimated, and in fact non-diffusive models have been proposed in the mathematical literature, see also [START_REF] Agutter | Intracellular transport mechanisms: a critique of diffusion theory[END_REF] for an interesting criticism of the idea of diffusion in intracellular transport mechanisms. Also the model should be made more realistic, for instance by the addition of an export mechanism for cargo and introducing a more complete cell environment. In our opinion, the purpose of this kind of mathematical models is as well that of proposing specific experimentation, not yet available in the current literature, such as single molecule tracking and localized cargo activation for estimating the needed parameters, which could greatly improve our understanding of the cell dynamics.

Figure 1 .

 1 Figure 1. Reaction scheme of Importin-mediated transport between cytoplasm and nucleus.

  Rt + Tc → Tr + C Ω n r3(Rt, Tc) := k3 Rt Tc k3 0.1 (μM s) -1[START_REF] Riddick | A systems analysis of importin-α-β mediated nuclear protein import[END_REF]

Figure 2 .

 2 Figure 2. Simple model of the MTs network as an aster originating from a center located in the proximity of the nuclear envelope. The cell environment schematized here has been used in the numerical experiments below, assuming a cell diameter of 60μm.

Figure 3 .

 3 Figure 3. Tumor suppressor Rb uniformly distributed in the cytoplasm at activation. Concentration of Rb (C + T c + D c ) in the cell: with active transport (above) and inhibited binding (below).

Figure 4 .

 4 Figure 4. Accumulation of the tumor suppressor Rb in the nucleus.

Figure 5 .Figure 6 .

 56 Figure 5. 3D-simulation of Rb nuclear import. Concentration of some of the molecular species involved at t = 100s: with active transport (above) and inhibited binding (below). Notice the localized concentration of the transport receptor and cargo complex T c in the proximity of the MTs aster when the MTs are active.

Figure 7 .

 7 Figure 7. Accumulation of the tumor suppressor Rb in the nucleus after peripheral activation in 2D.

Figure 8 .

 8 Figure 8. Three-dimensional plot of the concentration in the cell (in 2D case) of T c (left) and D c (right) after 20 seconds from activation, corresponding to the second plot above in Figure 6

Figure 9 .

 9 Figure 9. Peripheral cargo activation. 3D-simulation of Rb nuclear import. Concentration of T c : with active transport (above) and inhibited binding (below).

Table

	RanGTP	R t	d r	22	[13]
	RanGDP (with NTF2)	R d	d r	20	[23, 56]
	cargo	C	d c	12	calculated
	receptor	T	d t	14	[13]
	RanGTP•receptor	T r	d tr	14	[13]
	cargo•receptor	T c	d tc	10	calculated

Ran model reactant variable diff. coef. value (μm 2 /s) reference

Acknowledgments

We want to thank Enrico Bersani for his work done in the preliminary part of this project. We also thank Patrizia Lavia and Giulia Guarguaglini (IBPM-CNR) for their patient and useful explanations, sometimes about basic principles of biology, and gratefully acknowledge Joseph Pawlowsky (Fullbright scholar at IBPM) for conducting preliminary experiments in support of the model in human cells.

Finally, this work was partially supported by the Interdipartimental CNR Project Bioinformatics, leaded by Luciano Milanesi (IBT-CNR).