key: cord-0003039-58xrrowk
authors: Knodel, Markus M.; Nägel, Arne; Reiter, Sebastian; Vogel, Andreas; Targett-Adams, Paul; McLauchlan, John; Herrmann, Eva; Wittum, Gabriel
title: Quantitative Analysis of Hepatitis C NS5A Viral Protein Dynamics on the ER Surface
date: 2018-01-08
journal: Viruses
DOI: 10.3390/v10010028
sha: 502fc7e74c22bd0763542416c09589b5dca35b75
doc_id: 3039
cord_uid: 58xrrowk

Exploring biophysical properties of virus-encoded components and their requirement for virus replication is an exciting new area of interdisciplinary virological research. To date, spatial resolution has only rarely been analyzed in computational/biophysical descriptions of virus replication dynamics. However, it is widely acknowledged that intracellular spatial dependence is a crucial component of virus life cycles. The hepatitis C virus-encoded NS5A protein is an endoplasmatic reticulum (ER)-anchored viral protein and an essential component of the virus replication machinery. Therefore, we simulate NS5A dynamics on realistic reconstructed, curved ER surfaces by means of surface partial differential equations (sPDE) upon unstructured grids. We match the in silico NS5A diffusion constant such that the NS5A sPDE simulation data reproduce experimental NS5A fluorescence recovery after photobleaching (FRAP) time series data. This parameter estimation yields the NS5A diffusion constant. Such parameters are needed for spatial models of HCV dynamics, which we are developing in parallel but remain qualitative at this stage. Thus, our present study likely provides the first quantitative biophysical description of the movement of a viral component. Our spatio-temporal resolved ansatz paves new ways for understanding intricate spatial-defined processes central to specfic aspects of virus life cycles.

When a virus hijacks host cells of the body, its aim is to subvert cellular metabolism in order to create new virus progeny, which will infect surrounding cells and disseminate the infection. models of the vRNA replication cycle of HCV. These models remained qualitative at this stage due to a lack of experimentally validated parameters [27] . Therefore, the present paper demonstrates the estimation of a parameter needed for spatio-temporal resolved diffusion-reaction models of the vRNA replication cycle of HCV. Even though we estimate only one parameter that is needed for such models, our study shows that the estimation of such parameters is a highly nontrivial task. This tasks asks for the combination of advanced experiments and simulations. Future studies are needed to estimate the diffusion constants of all agents entering models of HCV replication. These studies will ask for further experimental studies which have to been combined with simulations in the manner we present here.

Hence, the detailed spatially-resolved understanding of the NS5A dynamics paves a way to advanced quantitative spatially-resolved understanding of HCV replication dynamics. Such a novel approach forms the basis of subsequent approaches extended initially to the other components of the HCV vRNA replication cycle and later, applied to other virus systems. Our approach may be used as first step of computational-based quantitative research in virology using advanced spatio-temporal resolved modeling and simulations techniques. To our best knowledge, the approach to apply highly precise numerical methods upon realistic geometric setups is a technique that has not yet been applied within the computational virology field.

The in silico description of the movement of NS5A is built on three basic components:

• Experimental data for the dynamics: FRAP time series data [34] recorded the intracellular dynamics of NS5A movement [10] at perinuclear zones. These data report the dynamics of NS5A under a variety of biological conditions. Namely we consider two different cell types explained in detail below. • Geometric setup: We use previously published confocal microscopic microscopy z-stack data [5] of cells labeled with ER markers which allow for reconstructions of realistic ER surfaces. These fine level data provide the geometric constraints for NS5A movement. • A model and corresponding simulations: Our previous model of NS5A dynamics [35] has not been adapted to biological data so far. In this study, we perform simulations using an extended version of the model and fit the simulation parameters in order to match the experimental data.

In the below, we describe these components in greater detail i.e. the FRAP time series experimental data, realistic ER geometry reconstructions, and we provide a short summary on NS5A movement properties that underpin the model we describe as basis for the sPDE simulations of NS5A dynamics (introduced by us recently [35] ).

FRAP experiments [32, 33] rely on the bleaching and refilling of a region within a cell corresponding to the intracellular location of a fluorescing component of interest, in this case, NS5A. The protein is distributed (and fluoresces) homogeneously at the beginning of the experiment. In a first step, a strong laser beam deletes the fluorescence in a defined region; the cylindrical FRAP region of interest (ROI) F . Subsequently, fluorescence recovers due to the influx of surrounding proteins from the surrounding unbleached region U . Measuring the intensity increase in F by a repetitive application of a second (soft) laser beam, allows to track the movement dynamics of the labelled NS5A proteins.

FRAP experiments were previously conducted, cf. [34] . The 20 time series published in this paper [34] were divided into two groups of in vitro cell types, namely NS5A/Alone and NS5A/OtherNSPs [34, 36, 37] , with 10 time series for each of the groups.

"NS5A/Alone" are Huh7 cells transfected with a DNA construct that encodes a NS5A-GFP fusion protein, but no other virus proteins are present. "NS5A/OtherNSPs" are Huh7 cells that contain a HCV genotype 2a subgenomic that encodes a NS5A-GFP fusion protein together with HCV non-structural proteins NS3, NS4A, NS4B, and NS5B. Thus; NS5A/OtherNSPs cells constitutively replicate HCV RNA but do not make virus since they don't produce the virus structural proteins (SPs). They do authentically replicate the virus RNA synthesis machinery though. In NS5A/Alone cells, NS5A-GFP is more mobile because it is likely not sequestered into membranous web regions to the same degree as it is in NS5A/OtherNSPs cells. NS5A/Alone cells allow the investigation of transport processes for less restricted NS5A-GFP. In contrast, NS5A/OtherNSPs cells have NS5A-GFP expressed as part of NS3-5B where NS5A diffusion is more restricted to replication sites.

NS5A anchors directly to the ER surface after its cleavage from the polyprotein [7, 28, 34] (like all HCV NSPs [1] ). Hence its movement is restricted to the ER surface and/or membranes derived from this organelle. As a first approximation, we may assume that the movement of this "free" NS5A can be modeled as diffusive movement. However, a proportion of NS5A can subsequently cluster to form vesicles, such as double membrane vesicles (DMVs) [6] . These vesicular structures contribute to the design of the membranous web. To our best knowledge, the relationship between the subpopulations of NS5A that cluster to form DMVs and NS5A that remains free is not yet understood. This applies equally to the NS5A/Alone data where only NS5A is expressed alone, and also for the NS5A that is expressed together with other NSPs, i.e., the NS5A/OtherNSPs cell type. Thus, the degree at which NS5A is associated with DMVs may be different for the two populations of NS5A we describe in the present study i.e., NS5A/Alone compared to NS5A/OtherNSPs.

At peripheral cell regions, the DMV clusters have a tendency to show random walk properties, i.e., to "jump around" [8] [9] [10] 38] . The peripheral NS5A foci (which are putative replication complexes) are frequently highly motile and capable of rapid long-range traffic. Hence, advanced models of NS5A will need to incorporate and addresses the motile populations of NS5A foci, for which particle tracking analysis will be an important evaluation tool, cf, e.g., [39] .

However, this "jumping" movement characteristics of NS5A clusters does not appear at perinuclear regions. Indeed, FRAP analysis of intense perinuclear NS5A foci have revealed a relatively static internal architecture [9, 34] . Since the FRAP experiments of NS5A [34] (which are the experimental basis for this study) were performed at perinuclear regions, we did not take into account the highly motile subpopulations of DMVs at peripheral subcellular locations during our modeling. Therefore, our present study focuses upon NS5A properties within perinuclear regions of the cell. Extension of the model to include peripheral NS5A properties [8] [9] [10] 38] will be the subject of forthcoming analyses and will require particle tracking analysis and single particle movement to complement our present diffusion/continuum model based framework.

First FRAP-based approaches were performed about 40 years ago by e.g., Soumpasis [40] and Axelrod et al. [41] . These analyses were based upon the assumption of pure 1D/2D processes, i.e., of mobility which takes place unrestricted within the 2D plane such that in many cases, the models can be reduced to 1D cases, i.e., exponential descriptions allow for solving the models. These basic assumptions are also a central part of more recent approaches for the modeling of FRAP experiments, cf. e.g., [42] [43] [44] [45] [46] . and also for similar modeling approaches for other fluorescence experiments like FLIP (fluorescence loss in photobleaching) [47] .

Our aim is to model FRAP experiments of the NS5A protein. Diffusive movement of NS5A is restricted to the ER surface while a portion of NS5A clusters to DMV complexes, which likely originate from the ER surface. The scope of this study is to reproduce the experimental time series from data previously published, cf. [34] . These data were produced within perinuclear regions where DMV complex traffic is not a common phenomenon.

We start deriving a transport model for NS5A on the ER surface for a given configuration E (i.e., ER geometry). In the model, we assume that it is reasonable to distinguish between two groups of NS5A molecules:

Let c (−) and c (+) denote the surface concentrations of NS5A occuring freely (−) (i.e., not clustering to DMVs or web regions) and in clustered form (+) (i.e., accumulated to DMVs respectively web regions) respectively. Note, that this definition tacitly assumes that it is reasonable to average over a certain area, which is large compared to the representative size of a cluster of molecules for either of the two types. Moreover, it is assumed, that transport of the two NS5A species is due to diffusive transport along the ER, i.e., in E the equations

hold. Here ∆ (T) denotes the Laplace-Beltrami-Operator, which is the projection of the Laplace operator to the tangential space of the two dimensional ER-hypersurface E that is embedded into R 3 [48] . The right hand side is a first order transformation model for mass transfer between both (±)-species. This model is further simplified as follows: (i) transport of the (−)-type is much faster than of the (+)-type, and (ii) the characteristic time of the reactions is much faster than the characteristic time for (diffusive) transport, i.e.,

Here, λ is the characteristic length for the experiment (i.e., the diameter of the blank spot generated by the laser beam in the FRAP experiment). Under these assumptions, we replace (1) by the quasi-steady-state equilibrium

and obtain a reduced version for (2):

Note that this can formally deduced by summing up (1) and (2) and approximating D + ≈ 0 due to (3) . Introducing the effective concentration

and the (observed) effective diffusion coefficient

yields the central sPDE for the dynamics of NS5A fluorescence expressed in terms of the intensity density

on the ER surface domain:

This is the "master equation" we use for the simulation of the spatially resolved dynamics of the data. In contrast to (5) , this equation includes the intensity reduction caused by the intensity measurement with the soft laser as a "pseudo reaction" by a constant r p . This effect is discussed in greater detail in the forthcoming Section 2.5. Note, that the constant r p does not correspond to a degradation of NS5A itself. This occurs at a much bigger time scale and is not included in the model.

Note that due to the presence of two different cell types, our framework models two different concentrations for each one of these cell types. Whenever it is important to distinguish, we indicate this by an addition index. We write c (±) A and i A ns5a for concentrations and intensities of cell type NS5A/Alone, and c (±) N and i N ns5a for the same quantities when referring to cell type NS5A/OtherNSPs. Likewise, we later derive diffusion constants D A ns5a , D N ns5a and measurement induced intensity reduction rates r A p , r N p for both cell types. The additional indices for the cell lines will be suppressed if the forthcoming statements are valid for both cell types.

In order to match the experimental setup, initial conditions for (8) are provided independently for bleached and unbleached zones/subdomains respectively. (The upper indices refer by now to the subdomains and not to the cell type, the cell type index is suppressed as explained before, since the explanations are valid for both cell types).

Due to the bleaching of the FRAP ROI with the strong laser at the beginning of the experiment, the intensity of the FRAP region is much smaller than of the surrounding unbleached region,

The initial values of the simulations, i F 0 and i U 0 , are determined by experimental values according to i F 0 = I F 0 /I 0 and i U 0 = I C 0 /I 0 . Here, I 0 denotes the intensity in the FRAP region before bleaching, I F (t 0 − ), where 0 < → 0. After bleaching, the intensities I F (t 0 ) and I U (t 0 ) are determined independently. The value for I U (t 0 ) is determined in an unbleached control region C ⊂ U . In this region, intensity is continuously monitored by independent application of a soft laser. The control region C is located in large distance from F , such that the mutual influence is negligible. All values are determined after subtraction of the background noise.

The sPDE (8) computations were performed with UG4 [49] [50] [51] based on Finite Volumes discretisations [52] [53] [54] and massively parallel multigrid solvers [55, 56] . (UG has been used successfully within various areas of computational physics [57, 58] and biophysics, namely computational neuroscience [59] [60] [61] [62] and computational pharmacology [63, 64] ). For technical details concerning the discretization methods and the massively parallel multigrid solvers, we refer to our former paper [27] .

The application of the soft laser induces an intensity decay, which is reflected by the constant r p in (8) . Its value is determined by the signal decay in the (unbleached) control ROI C, assuming the following ODE:

Here, i 0 and r p were determined for each one of the 20 TMS and averaged for both cell types separately by fitting to the corresponding experimental values [34] . Figure 1 provides an example of the solution of (10) fitted to the experimental values of the control ROI for one particular time series. Aggregating these results for both cell types yields the distributions shown in Figure 2 . (The result for one TMS (NS5A/Alone #2) yielding a negative r p was discarded due to an assumed measurement error. Note: Counting of time series starts from 0 ranging to 9, hence e.g., the 10th TMS has index # 9). The stochastic distributions and averages were computed with standard algorithms [65] using R [66] .

The data shows that the values for r p are not correlated and can be distinguished. The reduction rate for NS5A/OtherNSPs cells is higher than for NS5A/Alone: 

Our aim was to simulate NS5A diffusion upon real cell geometries. 20 NS5A FRAP time series from a former study [34] are the experimental basis of this study. These 20 time series reflect the NS5A dynamics. Each time series corresponds to one experiment. Each experiment was performed with one in vitro cell. This means that each time series corresponds with one special cell geometry. Each time series corresponds to one cell structure. We wanted to solve sPDE (8) upon realistic reconstructed ER surfaces. The computational domain E for sPDE (8) corresponds to the ER surfaces for each one of the 20 time series [34] . We wanted to reconstruct the cell structure of each time series. On this way, 20 ER geometries would have arisen.

However, there is a general technical limitation: An in vitro cell can be used either for FRAP experiments, or for staining purposes. Simultaneous staining of compartments of a cell and recording of FRAP time series at the same cell is not possible. This means, one can use the cell either to record time series. Or one can use the cell to stain special compartments like the ER. Such staining is the basis for reconstructions. However, a cell used to perform a FRAP time series cannot be stained. Compartments which are not stained cannot be reconstructed. Since it is not possible to use the same cell to record dynamics and to stain compartments, the cell geometry of a FRAP time series cell cannot be reconstructed. To this end, the ER geometry of a FRAP time series cannot be used as basis for FRAP simulations. There was no possibility to reconstruct the ER geometries which correspond to the cells of the 20 FRAP experiments. There was no possibility to evaluate sPDE (8) upon one of the cell geometries where the experimental time series [34] were recorded.

Despite all these adversities, our aim was to perform our simulations upon realistic geometric environments. Therefore, we reconstructed 5 realistic ER surfaces [35] using NeuRA2.3 [67, 68] based on (with Huygens [69] deblured) z-stacks from another former study [5] . For each ER geometry, we were probing for 2 different exemplary FRAP ROIs (selected manually) with a size of 38 (µm) 2 as in the FRAP experiments [34] . Figure 3 shows an example of such a geometric setup. We performed various tests to ensure that the choice of these geometries is justified, cf. Sections 2.7 and 3.3. The basic geometry details (like vertex numbers of the unstructured triangular grids-about 10 6 at base level), discretization and solver properties were presented in our former paper [35] as basis for extended refinement stability checks of single sPDE evaluations (using heuristic values for the diffusion constant). The Appendices B and C give additional information on the geometries. (ER geometry indexing starts from 1 ranging to 10).

We had to estimate the NS5A diffusion constant for two different experimental scenarios, namely NS5A/Alone cells and NS5A/OtherNSPs cells, cf. Section 2.2. For the comparison with experimental FRAP ROI intensity (cf. Section 2.4), we define the integrated normalized luminosity

obtained from the simulations. i ns5a is computed as solution of the sPDE (8) for each one of the following combinations: Combining 10 reconstructed geometric setups with 20 experimental time series, yields a total of 200 combinations. In order to find optimal values of D ns5a for each one of the 200 combinations of time series-geometric setup, we used the Gauss-Newton procedure [70] .

The averaging for final diffusion constant values was performed over the single estimated values for the combination of all geometric setups for each cell type over its 10 respective time series. All averaging procedures and all distribution computations elaborated in this study were performed with standard algorithms [65] implemented within the program code R [66] .

Our diffusion model of NS5A on the ER surface as derived in Section 2.4 was applied to the experimental data based reconstructed ER geometries which we explained in Section 2.6. A screenshot of a single FRAP experiment simulation (8) for one of the combinations (time series-geometric setup) is shown in Figure 4 . Figure 4 is a screenshot from supplemental movie "S1 Video in supplementary material". Figures 5 and 6 depict typical comparison curves of experiment and simulation, i.e., of FRAP ROI intensities (as explained in Section 2.7). Throughout this paper, we use for figures the notations: "TMS"-"time series", "geo(m)"-geometry, "psr"-"pseudo reaction".

The above-mentioned sPDE techniques were applied to estimate the NS5A diffusion constant on the ER surface for the experimental time series of both NS5A/Alone and NS5A/OtherNSPs Huh7 cells.

In detail, the single simulations were performed as explained before in Section 3.1. The simulation intensity of the FRAP region (12) of these single in silico processes was then fitted to the experimental time series to obtain the optimal diffusion constant for each single combination of time series and ER geometry as explained in Section 2.7, for examples, cf. Figures 5 and 6 .

The final diffusion constant was calculated as the average of the estimated diffusion constants of all combinations time series-geometric setups (with the methods described in Section 2.7). Before reporting the entire final values, we want to explain the investigations we performed to ensure the validity of our results.

To minimize the risk of artificial errors caused by the independent geometric setups, we tested intensively the influence of the time series and of the geometric setups on the final averaged diffusion constant. We found only a minor influence of the geometric setups (for both cell types). Figure 7 shows the distributions of the single and averaged values. The thick distribution curves of Figure 7a reflect the results for the single geometries as averaged over all time series, whereas the thick distribution curves of Figure 7b reflect the results for single time series averaged over all geometries. The distribution curves show that the geometric setups have only a minor influence in comparison to the time series. Therefore, the use of other realistic geometric setups from independent ER staining data [5] instead of the non available original ones of the TMS data [34] seems to be justified. 11)) which models the measurement process induced signal reduction within (8) . The in silico curves are that curves which arise for the estimated optimal value of D N ns5a for the time series # 5 of the NS5A/OtherNSPs cell case adapted to the reconstructed ER geometry E 1 .

The numerical stability of the results averaged over the geometric setups and the time series for the respective two cell types concerning the refinement level was intensively tested and showed sufficiently converging results for one-fold spatial refinement based sPDE evaluations. For details, cf. Appendix E.

The diffusion constant estimations were also done for zero pseudo reaction demonstrating the importance of the consideration of the measurement based intensity reduction. Only the NS5A/Alone cell case allowed for the derivation of reliable results for r p = 0. For the detailed results, we refer to Appendix D. Numerically, we also tested for the dependence of D ns5a on r p (covering also non-biophysical values) and found an excellent linear agreement, cf. Appendix F.

Comparative computations were done for a continuum model based on a planar 2D geometry with circular F to test for the influence of the curved ER manifold on the dynamics. This simulation type corresponds to the classical type of FRAP modeling, because the geometric structure is not resolved in this simplified case. The fit for each one of the 20 time series and the averaging for NS5A/Alone and NS5A/OtherNSPs cells showed a significant difference compared to the case when the ER structure is resolved. This demonstrates the importance of the use of a correct model [71] . Figure 8 depicts a screenshot of a simulation at the 2D plane geometry referring to supplemental movie "S3 Video in supplementary material". In Appendix C we show the distribution of the estimated diffusion constant corresponding to the single time series and the averaged results of both cell types (in a similar way as performed in Figure 7 for the ER surface case). (For example, the left y value 2.5 corresponds to ER geometry E 3 and TMS # 5). Half thick symbols (shown in the middle of each ER geometry region) correspond to the averaged values over all TMS D ns5a | T for the respective ER geometry as reported in Section 2.6. Aggregating these averages over all ER geometries yields distributions (thick continuous lines, scale shown on right y axis). Thick symbols (shown on top) correspond to the averaged values D ns5a . Note: The total averages are identical for (a,b) and are reported in Table 1 . (Red high concentration, blue low concentration). At the beginning, the bleached region has low concentration, but the 2D diffusion refills it again during the process, i.e. the color shifts slowly to red again because from the high concentration unbleached region around, fluorescating NS5A is diffusing inside. The ER structure is neglected.

The detailed investigations we performed and which we describe previously demonstrated: We have derived stable averaged results of the optimized NS5A diffusion constant on the ER surface for the NS5A/Alone and NS5A/OtherNSPs cell types based on 10 realistic reconstructed ER scenarios combined with 10 respective experimental FRAP time series.

The final averaged results for both cell types are shown in Table 1 for the ER manifold surface and for the 2D planar case.

The diffusion constant for the NS5A/Alone cells was approximately 4-fold larger than that of the NS5A/OtherNSPs cells (for both geometry types, i.e., ER surface and 2D planar geometry case).

Hence, the estimated value for the NS5A/OtherNSPs is significantly smaller then that one of NS5A/Alone and indicates that in the presence of other NSPs, the NS5A mobility is substantially reduced, presumably due to a higher amount of NS5A clustered to DMVs or membranous web regions.

The parameter estimations based on simplified 2D planar geometry, rather than the ER surface setup, caused a decrease of the diffusion constant values by a factor of approximately 2 (for both cell types). Thus, geometric simplifications as used often within simulations of biophysical processes change the results significantly. Table 1 . Averaged final NS5A diffusion constant D ns5a as described in Section 3.7. The final values are computed by means of the averaging process of the single results (as described in Section 2.7) which was shown graphically within Figure 7 . We give results for the case of the use of the ER geometry setups as described in Sections 3.1 and 3.2, but also for a simplified classical 2D planar consideration, cf. Section 3.6. 

We derived values for the NS5A diffusion constant of two experimentally important cell types, namely NS5A/Alone and NS5A/OtherNSPs cells.

The results show a 4-fold decrease in the diffusion constant observed within NS5A/OtherNSPs compared to NS5A/Alone cells (cf. Table 1 ). This decrease in the diffusion constant indicates a substantial effective reduction of the mobility of NS5A when the protein is expressed together with other NSPs in comparison to the situation when NS5A is expressed alone. Since the relation between those NS5A proteins clustering together to form DMVs or membranous web regions and those which remain freely diffusing on the ER surface has not been quantified experimentally to date, we assume that the diffusion constant of the freely diffusing NS5A species is identical for both cell lines. Therefore, the likely reason for the decreased effective diffusion constant observed when other NSPs are present is the presence of other NSPs enhancing the accumulation/clustering rate of NS5A into DMVs or membranous web regions. The clustered NS5A protein species have greatly reduced mobility. If more NS5A proteins cluster, the overall effective diffusion constant is reduced. This occurs when other NSPs are expressed (as in the case of our NS5A/OtherNSPs cell type) and not when only NS5A is expressed alone (which is valid for our NS5A/Alone cell type).

We want to analyze this property quantitatively based upon our mathematical model: According to the model expressed in Section 2.4, D − identifies the diffusion coefficient of the non-clustered (−) fraction of NS5A. Since this species is freely diffusing along the ER surface, we can assume that this constant is identical in both cell types, i.e., D

N . The four-fold increase of the effective diffusion constant once NS5A is expressed alone indicates obviously that less NS5A belongs to the clustered (+) fraction of NS5A. Even though some fraction of NS5A clusters also for the NS5A/Alone case, this fraction is substantially reduced in comparison to the NS5A/OtherNSPs case. Therefore, even though c (+)

N . Therefore, we may consider the limiting case c (+) A −→ 0 as a first order approximation for mathematical-quantitative consideration. Hence, we want to consider the most extreme case assuming for a moment that all NS5A do not cluster if NS5A is expressed alone. Of course, this is likely not the actual case, but the strongly enhanced portion of NS5A which is not clustering in the NS5A/Alone cell case compared to the NS5A/OtherNSPs case allows to consider this extreme approximation as some sort of first order approximation. The first-order hypothesis that all NS5A are not clustering for the NS5A/Alone case indicates that k A on −→ 0 within (2). Since we assume that the diffusion constant of both species of NS5A is the same for both cell types, D

N , the difference for the both cell types originates within the different binding rates, k A on = k N on respectively k A off = k N off . Within our modeling framework, this hypothesis corresponds to differences in D ns5a from Equation (7): Under the assumptions (3), the fitted results to ratio k N on /k N off = 3 for long times for the limit case k A on 0. That indicates a ratio of 1:3 for c (−) N to c (+) N due to (4), i.e., for free versus clustered species for the NS5A/OtherNSPs case within our first order approximation. This means that in the NS5A/OtherNSPs case, an appreciable amount of NS5A molecules are clustered in membranous web regions / DMVs, and thus relatively immobile.

The diffusion of NS5A species on the ER surface that are not clustered into DMVs/membranous webs are likely not completely "free" on the ER surface due to interactions with a huge possible number of host interacting proteins (about 130 [72] ). However, our results indicate that the most important obstruction of the NS5A movement originate from its clustering to DMVs or web regions when additional NSPs are present. This would account for the 4-fold diffusion constant decrease once other NSPs are expressed.

To the best of our knowledge, applying spatio-temporal modeling to biological data has not been previously performed in this manner and forms basis for further detailed investigation and refinement using subsequent models [27] .

We are developing spatio-temporal resolved models of the HCV replication cycle within single liver-derived cells. In our previously published paper [27] , we developed spatially resolved models of the vRNA cycle. These models considered the major components of the vRNA cycle, namely vRNA, NSPs and a host factor. The model was qualitative rather then quantitative due to a lack of experimentally-derived parameters. The diffusion constant for the NS5A/Alone case which we derived in this study is a candidate for the class of diffusion-reaction models of HCV replication [27] . Future work will ask for the determination of the diffusion constants of the other components of the vRNA replication cycle.

The control computations using a 2D continuum model gave values which for the two cell types were about a factor 2 smaller than for the ER surface simulations. The differences between ER manifold sPDE and 2D continuum PDE results are in agreement with observations based on particle based PDE evaluations [73, 74] . Only few publications deal with the quantitative analysis of spatio-temporal properties of virus proteins [31, 75] at all. Firm values for viral protein diffusion constants have not been reported in the literature. For metabolism proteins and other physiologic ingredients, some evaluations exist, cf. e.g., [73, 74, [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] [86] [87] . Besides the studies published within [73, 74] , these approaches do not take into account the detailed structure of the ER. In most cases, the PDE evaluations are based on simplified techniques. Advanced numerical methods have been applied only in very few cases within the field of cellular simulations at the ER [88] yet.

The estimation of biophysically meaningful results for the diffusion constant of the important NS5A viral protein on the curved ER surface manifold is an intellectually-stimulating contribution to the young field of spatio-temporal resolved research within computational virology. For the first time, the derived results give a quantitative biophysical description of the movement properties of a crucial viral protein related to various steps of virus replication-indeed, we believe the parameter estimated by us represents the first quantitative description of the movement characteristics of any viral component at an (intra)cellular level. These results are intended to enter spatio-temporal resolved biophysical models of HCV replication in the type we have presented previously in our model paper [27] , and, later, the techniques can be applied to similar systems [2] [3] [4] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] 89] including non-viral processes [88] , e.g., metabolic protein [73, 74, [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] [86] [87] dynamics.

Our novel approach to introduce spatio-temporal resolved simulation techniques into computational virology paves a way for previously unexplored detailed biophysical understanding of virus replication dynamics, for example to unveil the relationship of form and function, as we have demonstrated already within our former paper [27] , or to reveal areas of the virus life cycle amenable to novel antiviral intervention that conventional biology may miss e.g., spatial dependence of virus-encoded factors within specfic intracellular regions. In the future, this avenue of research has the possibility of substantially impacting our understanding of complex biological systems.

Supplementary Materials: The following are available online at www.mdpi.com/1999-4915/10/1/28/s1. Supplemental Movie S1 Video: Movie of FRAP simulation at ER geometry I, Supplemental Movie S2 Video: Movie of FRAP simulation at ER geometry IV, Supplemental Movie S3 Video: Movie of classical FRAP simulation at 2D continuum plane. (cf. also the Appendices A.1-A.3 for brief descriptions of the simulation movies). SVI, Netherlands) for his very friendly support in Huygens usage, backgrounds, and licensing. The HLRS Stuttgart is acknowledged for the supplied computing time on the Hermit and Hornet super computers [90] , and Michael Lampe for very friendly technical support on the G-CSC cesari cluster. The authors acknowledge the Goethe Universität Frankfurt for general support and computational resources and the Politecnico di Torino for general support. This work has been supported in part by the "Fondazione Cassa di Risparmio di Torino" (Italy), through the "La Ricerca dei Talenti" (HR Excellence in Research) programme. The Authors wish to express their sincere thanks to the anonymous Referees for their thorough and critical reviews of our work. 

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

The following abbreviations are used in this manuscript: Simulation of NS5A FRAP experiment simulation at 2D continuum plane (computed by means of PDE description of NS5A diffusion with UG4, ∆ (T) → ∆, i.e., "trivial" 2D diffusion instead of manifold surface diffusion, use of "common" 2D Laplace operator ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 in (8)). Figure A1 depicts the basic ER geometries E i , i = 2, 3, 4, 5, i.e., all ER geeometries with one FRAP ROI example (ER geometry E 1 is shown in Figure 3 ). Figure A1 . The reconstructed ER geometries E i , i = 2, 3, 4, 5 and exemplar FRAP regions. FRAP ROIs indicated in red, covering a surface of about 38 µm 2 enabling the reader to estimate the total size of the ERs. Table A1 depicts the number of the degrees of freedom (DoF) and the number of the basic elements (triangles) of the ER geometries corresponding to the different subdomains. Roughly, the number increases about a factor 4 for each refinement level. For more details, cf. [35] .

ER geometries 6-10 harbor the same number of DoFs and faces, however slightly different numbers within F / U caused by the shift of the FRAP ROI. Table A2 depicts the number of DoFs for the simple planar 2D geometry based classical 2D FRAP modeling which is used for comparison computations, i.e., when no ER structure is taken into account. cf. Section 3.6. Figure A2 depicts the distribution of the estimated diffusion constant corresponding to the single time series and the averaged results of both cell types.

For the NS5A/Alone cell case, we estimated the diffusion constant also setting r A p = 0 within the basic sPDE (8) . (For the NS5A/OtherNSPs case, this was not possible reliably, i.e., using r N p = 0 did not allow to estimate reliable diffusion constants). Figure A2 inherits the single (i.e., for each time series) estimated diffusion constant values, their distribution and the averaged diffusion constants evaluated for the 2D case. Figure A3 shows the single and averaged results and the distributions for the computations on the ER surface. Table A3 depicts the final averaged values for the case r A p = 0 for the NS5A/Alone case, for the ER surface based computations and for the classical 2D plane based computations. Table A1 . DoF number and total face number of the ER geometries at level L, extracted from our former paper [35] . 

In our former paper [35] , we investigated the refinement stability of single sPDE computations for heuristic diffusion constant values (setting r p = 0). We tested the stability of the temporal evolution of the FRAP ROI integrals I(t), (12) , under variation of spatial and temporal refinement. Therefore, we evaluated the relative differences under spatial level resp. time step size refinement.

Here, we analyze the refinement stability of the estimated averaged diffusion constant. As shown in our paper [35] , the influence of the refinement stability increases when the diffusion constant decreases. Since the pseudo reaction does not depend on spatial refinement, we only needed to check for the case r p = 0 in (8) in principle. Therefore, we check the NS5A/Alone cell case without pseudo reaction (r A p = 0). (The corresponding single and averaged results were discussed within Appendix D afore). Since the NS5A/OtherNSPs cell case did not allow for a parameter estimation without incorporating the pseudo reaction, we perform the check in this case with pseudo reaction (r N p = 0)).

In order to save computation time, we omitted to perform the level 2 computations for the second location of the FRAP ROIs and compared the results for one FRAP ROI per reconstructed ER geometry. Table A4 presents the investigations of the spatial refinement stability of the averaged Gauss-Newton results for the NS5A/Alone cell case and for the NS5A/OtherNSPs cell case (using respectively only one FRAP ROI per geometry for simplicity). Our standard time step size 0.1 s has only minor influence on the numeric error compared to the spatial refinement influence [35] . Hence, there was no need to check the stability of the averaged results concerning the time step size refinement. Finally we omitted extended refinement stability checks for NS5A/OtherNSPs cells combined with the planar 2D geometry since no additional insight would have been gained.

The averaged values considering only 5 geometric setups (Table A4 ) are negligible different (3%) from the final values averaged over all 10 geometric setups, Table 1 (for level 1 computations).

Based on the comparisons shown in Table A4 , we conclude that level 1 is sufficient for the ER computations and level 8 for the 2D case computations. Figure A2 . Classical FRAP analysis based upon the 2D computations explained in Section 3.6: Averages for diffusion constant estimation D ns5a -analyzed separately for NS5A/Alone (red-D A ns5a ), NS5A/OtherNSPs (green-D N ns5a ) as well as for NS5A/Alone setting r p = 0 (blue-D A ns5a | r p =0 ), i.e., assuming that the measurement process induces no intensity reduction. Thin points and error bars correspond to estimated values for D ns5a (shown on x axis) for a single TMS (indicated on the left y axis). Aggregating over all TMS yields distributions (continuous lines, scale shown on right y axis). Thick symbols (shown on top) correspond to the averaged values D ns5a reported in Table 1 for the 2D planar case and in Table A3 for the 2D planar case using r A p = 0 for the NS5A/Alone cell case. (for the NS5A/Alone cell case) as described in Section 3.7, but assuming vanishing signal reduction due to the measurement process, i.e., using r A p = 0 in (8) . The final values are computed by means of the averaging process of the single results (as described in Section 2.7) which is shown graphically within Figure A3 . We give results for the case of the use of the ER geometry setups as described in Sections 3.1 and 3.2, but also for a simplified classical 2D planar consideration, cf. Section 3.6. Values only given for the NS5A/Alone case, since reliable estimations were not possible for the NS5A/OtherNSPs case once we used r N p = 0. Figure A3 . Averages for NS5A diffusion constant D A ns5a estimation on the ER surface as described in Sections 2.2 and 3.3 for the NS5A/Alone cell case, neglecting the measurement process induced signal reduction, i.e., setting r A p = 0 in (8) . Simular analysis as in Figure 7 . Thin points and error bars correspond for (a,b) to estimated values for D ns5a (shown on x axis) for the combination of single TMS with single geometries (indicated on the left y axis, note different combinations for geometry and TMS). Aggregating over all TMS and ER geometries yields distributions (thin continuous lines, scale shown on right y axis) which are identical in both cases. (a) Each "row" corresponds to one time series combined with all ER geometries. (For example, the left y axis value 2.5 corresponds to the combination of TMS # 2 and ER geometry E 6 ). Half thick symbols (shown in the middle of each time series region) correspond to the averaged values over all geometries D ns5a | G for the respective TMS. Aggregating these averages over all TMS yields distributions (thick continuous lines, scale shown on right y axis). Thick symbols (shown on top) correspond to the averaged values D ns5a . (a) Each "row" corresponds to one ER geometry combined with all TMS. (For example, the left y value 2.5 corresponds to ER geometry E 3 and TMS # 5). Half thick symbols (shown in the middle of each ER geometry region) correspond to the averaged values over all TMS D ns5a | T for the respective ER geometry as reported in Section 2.6. Aggregating these averages over all ER geometries yields distributions (thick continuous lines, scale shown on right y axis). Thick symbols (shown on top) correspond to the averaged values D A ns5a . Note: The total averages are identical for (a,b) and are reported in Table A3 . Table A4 . Refinement stability: Averaged NS5A diffusion constant for the simplified planar 2D case and the ER surface manifold computations (using only one FRAP ROI per ER geometry. Nota bene: the afore reported final results cover all geometric setups, i.e., the final results reported in Table 1 are based upon the use of two FRAP ROIs per reconstructed ER geometry). Evaluation for different spatial refinement levels R of base geometry and relative change C in comparison to level before in percentage. Values for the NS5A/Alone cell case without pseudo reaction, i.e., r A p = 0 in (8) and for the NS5A/OtherNSPs cell case with pseudo reaction, r N p = 0. 

We tested for the influence of the numerical variation of the pseudo reaction constant on the final result of the diffusion constant, i.e., we analyze how the final diffusion constant changes under variation of the value of r p in our "master sPDE" (8) .

Therefore, we estimated the diffusion constant D = D ns5a also for (non-biophysical) values of r p in the region of the biophysical values. Hence, the estimations of D ns5a were performed not only for the estimated values of r p , but also for other ones close to them to test for the dependency.

Due to the smallness of the scale of r p -O 10 −3 -the variation of D can be very well approximated by a linear fit, since higher terms in the Taylor approximation are negligible. Using the linear regression ansatz D(r p ) = D 0 + f r p + O(r 2 p ) (A1)

we derived the factors D 0 , f . A graphical representation of the fit for the ER surface computations is given in Figure A4 . Figure A5 depicts the variation for the 2D case in graphical representation. Table A5 shows the numerical coefficients of (A1) for the ER surface calculations. Table A6 depicts the coefficients for the 2D case. The graphs, coefficients of the fits and the fit procedure are shown here in the appendices since this investigation is not of biophysical interest, but moreover serves as an additional numerical stability check. The quality of the linear fits, i.e., the excellent linear agreement, can be considered as a further hint of the stability of our results.

Whereas Figure 4 depicts a screenshot of the simulations at ER geometry I, Figure A6 shows an additional simulation screenshot of the ER computations, namely at ER geometry IV. The corresponding movie is attached as supplemental movie "S2 Video in supplementary material". Figure A4 . Dependency diffusion constant of NS5A on ER surface on pseudo reaction constant, linear fit of (A1) for numerical stability test reasons. "estim": estimated values, "fit": linear fit. 

Replication of hepatitis C virus

Architecture and biogenesis of plus-strand RNA virus replication factories

Dengue virus and Hepatitis C virus-induced replication and assembly compartments: The enemy inside-Caught in the web

Composition and three-dimensional architecture of the dengue virus replication and assembly sites

Visualization of double-stranded RNA in cells supporting hepatitis C virus RNA replication

Three-dimensional architecture and biogenesis of membrane structures associated with hepatitis C virus replication

Small molecules targeting hepatitis C virus-encoded NS5A cause subcellular redistribution of their target: Insights into compound modes of action

Daclatasvir inhibits hepatitis C virus NS5A motility and hyper-accumulation of phosphoinositides

A dynamic view of hepatitis C virus replication complexes

Dynamic imaging of the hepatitis C virus NS5A protein during a productive infection

Oligomerization of Ebola Virus VP40 Is Essential for Particle Morphogenesis and Regulation of Viral Transcription

Inclusion Bodies Are a Site of Ebolavirus Replication

Assembly and budding of Ebolavirus

Membrane binding and bending in Ebola VP40 assembly and egress

Full-length Ebola glycoprotein accumulates in the endoplasmic reticulum

Structural dissection of Ebola virus and its assembly determinants using cryo-electron tomography

Severe Acute Respiratory Syndrome Coronavirus ORF6 Antagonizes STAT1 Function by Sequestering Nuclear Import Factors on the Rough Endoplasmic Reticulum/Golgi Membrane

Intracellular Events and Cell Fate in Filovirus Infection

Shedding Light on Filovirus Infection with High-Content Imaging

Live-cell imaging of Marburg virus-infected cells uncovers actin-dependent transport of nucleocapsids over long distances

The ER in 3D: A multifunctional dynamic membrane network

A perspective on modelling hepatitis C virus infection

Mathematical Modeling of Subgenomic Hepatitis C Virus Replication in Huh-7 Cells

Modeling Subgenomic Hepatitis C Virus RNA Kinetics during Treatment with Alpha Interferon

A Multi-Variant, Viral Dynamic Model of Genotype 1 HCV to Assess the in vivo Evolution of Protease-Inhibitor Resistant Variants

Replication vesicles are load-and choke-points in the hepatitis C virus lifecycle

3D spatially resolved models of the intracellular dynamics of the Hepatitis C genome replication cycle

Small molecule inhibitors of the hepatitis C virus-encoded NS5A protein

Essential Role of Domain III of Nonstructural Protein 5A for Hepatitis C Virus Infectious Particle Assembly

Review of boceprevir and telaprevir for the treatment of chronic hepatitis C

Role for TBC1D20 and Rab1 in Hepatitis C Virus Replication via Interaction with Lipid Droplet-Bound Nonstructural Protein 5A

From fixed to FRAP: Measuring protein mobility and activity in living cells

Advanced Fluorescence Microscopy Techniques-FRAP, FLIP, FLAP, FRET and FLIM

Mobility analysis of an NS5A-GFP fusion protein in cells actively replicating hepatitis C virus subgenomic RNA

Multigrid analysis of spatially resolved hepatitis C virus protein simulations

Enhancement of hepatitis C virus RNA replication by cell culture-adaptive mutations

Mutations in hepatitis C virus RNAs conferring cell culture adaptation

HCV RNA traffic and association with NS5A in living cells

Diatrack particle tracking software: Review of applications and performance evaluation

Theoretical analysis of fluorescence photobleaching recovery experiments

Mobility measurement by analysis of fluorescence photobleaching recovery kinetics

Kinetic modelling approaches to in vivo imaging

FRAP analysis of binding: Proper and fitting

Analysis of binding reactions by fluorescence recovery after photobleaching

Fluorescence recovery after photobleaching: Application to nuclear proteins

Identification of biomolecule mass transport and binding rate parameters in living cells by inverse modeling. Theor

Computational modeling of fluorescence loss in photobleaching

Differential Geometry: Curves-Surfaces-Manifolds

Software Framework ug4: Parallel Multigrid on the Hermit Supercomputer

A massively parallel geometric multigrid solver on hierarchically distributed grids

UG 4: A novel flexible software system for simulating PDE based models on high performance computers

Elliptic Differential Equations: Theory and Numerical Treatment

Finite-Volumen-Und Mehrgitter-Verfahren Für Elliptische Randwertprobleme

A generalization of the vertex-centered finite volume scheme to arbitrary high order

Multi-Grid Methods and Applications

Iterative Solution of Large Sparse Systems of Equations

Models and simulations of variable-density flow in fractured porous media

Simulation of salinity-and thermohaline-driven flow in fractured porous media

Synaptic Activity Induces Dramatic Changes in the Geometry of the Cell Nucleus: Interplay Between Nuclear Structure, Histone H3 Phosphorylation, and Nuclear Calcium Signaling

Synaptic bouton properties are tuned to best fit the prevailing firing pattern

A Three-Dimensional Mathematical Model of Active Signal Processing in Axons

1D-3D Hybrid Modelling-From Multi-Compartment Models to Full Resolution Models in Space and Time

In-silico model of skin penetration based on experimentally determined input parameters. Part II: Mathematical modelling of in-vitro diffusion experiments. Identification of critical input parameters

A comparison of two-and three-dimensional models for the simulation of the permeability of human stratum corneum

A Handbook of Statistical Analyses Using R

Nonlinear anisotropic diffusion filtering of three-dimensional image data from 2-photon microscopy

Inertia Based Filtering of High Resolution Images Using a GPU Cluster

Scientific Volume Imaging B.V., Huygens Comute Engine, Software

FRAP and kinetic modeling in the analysis of nuclear protein dynamics: What do we really know?

Serine phosphorylation of the hepatitis C virus NS5A protein controls the establishment of replication complexes

Effects of Organelle Shape on Fluorescence Recovery after Photobleaching

Simulations of (An)Isotropic Diffusion on Curved Biological Surfaces

Quantitative live-cell analysis of microtubule-uncoupled cargo-protein sorting in the ER

Dynamics and retention of misfolded proteins in native ER membranes

FRAP beam-size analysis to measure palmitoylation-dependent membrane association dynamics and microdomain partitioning of Ras proteins

Protein Diffusion in Mammalian Cell Cytoplasm

Dissecting the Contribution of Diffusion and Interactions to the Mobility of Nuclear Proteins

Transient Anomalous Diffusion of Telomeres in the Nucleus of Mammalian Cells

In Vivo Anomalous Diffusion andWeak Ergodicity Breaking of Lipid Granules

Dynamics of Single mRNPs in Nuclei of Living Cells

Nuclear body movement is determined by chromatin accessibility and dynamics

Macromolecular mobility inside the cell nucleus

Cajal Body dynamics and association with chromatin are ATP-dependent

Sensing Viruses by Mechanical Tension of DNA in Responsive Hydrogels

Single-Image Diffusion Coefficient Measurements of Proteins in Free Solution

Reaction diffusion modeling of calcium dynamics with realistic ER geometry

Serotype Specfic Differences in Dengue Virus Non-Structural Protein 5 Nuclear Localization

On Estimation of a Viral Protein Diffusion Constant on the Curved Intracellular ER Surface

We thank Konstantinos Xylouris (G-CSC) for very fruitful discussions on the evaluation of the simulation results, Martin Rupp (G-CSC) for very friendly technical support with respect to the implementation of the solvers, Ranjita Dutta-Roy (Karolinska Institute, Stockholm, Sweden) for profound explanations of FRAP experimental setup and data analysis, Stefan Groote (University of Tartu, Estonia) for proof reading the manuscript and very helpful hints, Tobias Denninger (Heidelberg, Germany) for very stimulating discussions about spatial resolution within computational biophysics, Alfio Grillo (Politecnico di Torino, Italy) and Jürgen Vollmer (MPI Göttingen, Germany) for very stimulating discussions about the subject, and Wouter van Beerendonk (Huygens