Keratin dynamics: modeling the interplay between turnover and transport.

Keratin are among the most abundant proteins in epithelial cells. Functions of the keratin network in cells are shaped by their dynamical organization. Using a collection of experimentally-driven mathematical models, different hypotheses for the turnover and transport of the keratin material in epithelial cells are tested. The interplay between turnover and transport and their effects on the keratin organization in cells are hence investigated by combining mathematical modeling and experimental data. Amongst the collection of mathematical models considered, a best model strongly supported by experimental data is identified. Fundamental to this approach is the fact that optimal parameter values associated with the best fit for each model are established. The best candidate among the best fits is characterized by the disassembly of the assembled keratin material in the perinuclear region and an active transport of the assembled keratin. Our study shows that an active transport of the assembled keratin is required to explain the experimentally observed keratin organization.

Introduction

The epithelial cytoskeleton is characterized by abundant keratin intermediate filaments (Fig 1). The cytoplasmic keratin filament network is responsible for the mechanical stress resistance of epithelial cells and contributes significantly to epithelial stiffness [1, 2]. The importance of keratins for epithelial tissue stability is reflected by a large group of genetic skin blistering diseases that are caused by point mutations in keratin-encoding genes [3, 4]. The mechanical functions not only rely on static resilience but also necessitate a high degree of plasticity, for example in migrating cells during wound healing [5]. The current view is that keratins act as general stress absorbers protecting epithelial cells not only against mechanical insults but also against irradiation or osmotic and microbial challenges. Thus, keratins are involved in heat shock response, apoptosis and organelle homeostasis [6]. Furthermore, functions affecting processes such as proliferation, differentiation and inflammation are also dependent on keratins (see recent reviews in [4, 7]).

Fig 1

Keratin network.

Image taken from a time-lapse fluorescence recording of a single hepatocellular carcinoma-derived PLC cell of clone PK18-5 [8] stably expressing fluorescent fusion protein HK18-YFP consisting of human keratin 18 and enhanced yellow fluorescent protein. Bar 10 μm.

All of these functions are tightly coupled to the keratin network dynamics (see S1 Video). Examination of cultured epithelial cells producing fluorescent keratins has provided evidence for the different mechanisms that are involved in the continuous renewal and reshaping of the keratin system [9]. On the basis of these observations, a biological model of the keratin cycle has been proposed by Leube and Windoffer et al. in [10, 11]. This biological model takes into account the assembly/disassembly and transport of keratins. In this study, it was proposed that assembly and disassembly occur in topologically defined regions with assembly taking place predominantly in the cell periphery while disassembly takes place primarily in the perinuclear region. The biological model further postulates active transport of insoluble assembly stages of keratins toward the nucleus and rapid diffusion of soluble subunits throughout the cytoplasm. To the best of our knowledge, how these different processes are coupled and regulated is not yet known.

In a previous work we developed methods to examine and quantify the keratin transport and turnover in epithelial cells [12]; the spatial distribution of the assembled keratin material in epithelial cells was available at 24 hours and 48 hours after seeding. The current effort is to use the available qualitative and quantitative observations to derive, from first principles, experimentally-driven mathematical models that could yield hypothetical predictions testable in laboratories. Our approach is unique, it translates experimental observations and data into a series of alternative plausible mathematical models or scenarios to further advance our understanding of the critical parameters in keratin cycling. Fundamental to this approach is the fact that optimal parameter values for each scenario are established and out of this set, a single scenario is identified that best fits experimental observations and data.

Hence, a collection of mathematical models resulting from different assumptions of the keratin transport, assembly and disassembly is designed to investigate the effects of the interplay between turnover and transport on the keratin organization. The collection is built as a well-designed scientific experiment by considering control and knockout of processes. To highlight and confirm the importance or existence of a given process, scenarios in which the process is absent are also considered and tested. Model responses are then compared to experimental observations and data published in [12] to identify optimal parameter values that yield the best fit of each of the models to the experimental data. Finally, we identify, using an information-theoretic approach, the best scenario or model given the data and candidate models under study.

By employing this reductionist phenomenological approach and systematic evaluation of the different scenarios we not only confirm the proposed transport features of the keratin cycle and the restricted disassembly in the perinuclear region but also find that the assembly throughout the cytoplasm fit best to the experimental data. Furthermore, our particular approach allows us to demonstrate that the inward motion experimentally observed is not an emergent behavior but is an inherent property of the keratin material organization and it is due to an active transport thereby confirming recent experimental observations [10, 11].

Methods

Experimental data

In Moch et al. [12], the spatial distribution of the assembled keratin material in epithelial cells is measured for 15 minutes at 24 hours and 48 hours after seeding. The shape of each epithelial cell is normalized to fit a circle of fixed radius [13]. The average normalized spatial distribution is calculated over 50 cells at 24 hours (resp. 84 cells at 48 hours). The average speed and direction of the motion of the assembled keratin material are measured and determined at every location within the normalized cell. Finally, at every spatial location, the net assembly/disassembly is calculated. Hence, regions with preferential assembly and disassembly are identified. We will refer to regions of assembly as Sources and regions of disassembly as Sinks. More details on the experimental data can be found in [12, 14].

In the present work, cells are represented as one-dimensional cross-section domains. A diameter of the normalized cell is used as the spatial domain which is centered at the center of the cell and is of length 2L with L = 22.5μm. Moch et al. recorded the fluorescence intensity of fluorescent protein-labeled keratins in cells. Assuming a proportionality between fluorescences and concentrations, and knowing from [15] the mean concentration for keratin in keratinocytes, we convert fluorescence intensities to concentrations (μM) as follows: Concentration = Fluorescence × (Mean Concentration / Mean Fluorescence). In Fig 2, the average spatial distribution, the speed of the assembled keratin, regions of assembly (Sources) and regions of disassembly (Sinks) on the one-dimensional cross-section domain are displayed.

Fig 2

Experimental data adapted from [12].

2.1: The spatial distribution of the assembled keratin material. The mean concentration of the keratin material used in the conversion of fluorescence intensities to concentrations is estimated to be equal to 520μM in [15]. Circles represent raw experimental data; curves are fits of experimental data. Details on P(x) and f (x) can be found in Appendix 4. 2.2: The speed of the assembled keratin material at 24 hours. 2.3: Regions of assembly and disassembly denoted Sources and Sinks respectively. In all figures, the cell is represented by a one-dimensional cross-section domain centered at the center of the cell; plasma membrane positions are at ±L = ±22.5μm, the nuclear envelope is located at ±7.5μm and the center of the cell is located at zero.

Mathematical models

To study its organization in cells, the keratin material is categorized into a soluble pool composed of the soluble keratin, and an insoluble pool representing the assembled keratin material. The state variables used to represent the soluble and insoluble pools are: From here onwards, the one-dimensional spatial domain representing the cell is defined by where x = 0 is the center of the cell and x = ±L are the boundary positions of the plasma membrane (Table 1). The general framework of the model, derived from first principles based on experimental observations, takes into account the turnover and transport for both soluble and insoluble pools and can be stated verbally as: Hence, the generalized model has the following expression; stated mathematically as: where T (⋅) (resp. T (⋅)) is the functional term describing the transport of the soluble pool (resp. insoluble pool). For example, (resp. ) where J (x, t) (resp. J (x, t)) describes the flux of the soluble pool (resp. insoluble pool) at position x from the left (x = −L) to the right (x = L) at time t. The function 𝓐(S) (resp. 𝓓(I)) is the assembly term (resp. disassembly term). To investigate the interplay between turnover and transport on the organization of the keratin material several assumptions are proposed for each of the functionals T , T , 𝓐 and 𝓓.

Table 1

Model parameters.

Parameter	Definition	Value (Unit)	Reference
L	Half-length of spatial domain Ω	22.5 (μm)	[12]
f 0(x)	Initial distribution of assembled keratin in cells at 24 hours	Eq (17) (μM)	[12]
D S	Diffusion coefficient of soluble pool	0.88±0.08 (μm 2/s)	[9]
D I	Diffusion coefficient of insoluble pool	≤ 10−3 D S (μm 2/s)
u	Speed of insoluble pool	0.002 − 0.008 (μm/s)	[12, 16–18]
v(x)	Almost constant speed of active transport of insoluble pool	Eq (12) (μm/s)	[12, 18]
u(x)	Variable speed of active transport of insoluble pool	Eq (13) (μm/s)	[12]
k ass	Rate of assembly of soluble pool	TBD (linear model: s −1, nonlinear model: μM.s −1)
k ass(x)	Localized rate of assembly of soluble pool	Eq (14) with k ass TBD (linear model: s −1, nonlinear model: μM.s −1)
k S	Concentration for half-saturation of assembly rate (nonlinear model)	TBD (μM)
k dis	Rate of disassembly of insoluble pool into soluble pool	TBD (linear model: s −1, nonlinear model: μM.s −1)
k dis(x)	Localized rate of disassembly of insoluble into soluble pool	Eq (15) or (16) with k dis TBD (linear model: s −1, nonlinear model: μM.s −1)
k I	Concentration for half-saturation for disassembly rate (nonlinear model)	TBD (μM)

The model parameters used in all the scenarios. (TBD = To Be Determined by fitting model solutions to experimental data).

S(x, t) denotes the concentration of the soluble keratin material at position x at time t,

I(x, t) denotes the concentration of the assembled keratin material at position x at time t.

Modes of transport

Molecules of the soluble pool are assumed to be subjected to the Brownian motion; the soluble pool is diffusible with a diffusion coefficient D . The functional term of transport for the soluble pool is given by Passive transport for both pools is assumed in all models. Diffusion is assumed for the insoluble pool to describe the wiggling motion of the keratin filaments in cells. The diffusion coefficient of the insoluble pool is set to be much smaller than that of the soluble pool: 0 < D ≤ 10−3 D . It is assumed that only the insoluble pool can be driven by an active transport. Experimental evidence show that the assembled intermediate filament proteins in the form of filament precursors, squiggles or filaments move along microtubules and actin filaments by interacting via molecular motors [11, 18–21]. An active transport (also called the inward drift) of the insoluble pool from the plasma membrane to the center of the cell is hypothesized based on reports of the motion of the assembled keratin material mostly towards the nucleus in epithelial cells [11, 17, 18]. The speed of the active transport is set to be almost constant v(x) everywhere or variable u(x) everywhere (Fig 3). Based on experimental observations, both speeds are assumed to decay around the nucleus towards the center of the cell. The variable speed u(x) is estimated using the profile of the average speeds measured in [12] (Fig 2.2). The magnitude of the almost constant speed v(x) is set to be the average value of the variable speed over the cell. Details on the derivation of the estimates of v(x) and u(x) are given in Appendix 1. Hence, the functional term for the transport of the insoluble pool can take three different forms: where v(x) is given in (12) and u(x) in (13); the two functions representing the speeds of the active transport are graphed in Fig 3. The function sgn(x) defined by describes the inward direction of the active transport at any location of the spatial domain Ω centered at zero.

Fig 3

Profiles of the active transport speeds u(x) and v(x) for the insoluble pool both of which are derived from the experimental measurements in [12].

The function u(x) is derived from the profile of speeds and v(x) is approximated as the average value of these speeds. Details on the derivation of u(x) and v(x) are given in Appendix 1.

Combining the modes of transport for the soluble and insoluble pools, three types of transport are considered for the keratin material in the epithelial cell.

Expressions of assembly / disassembly reactions

In the present work, the turnover is composed of two reactions: the assembly / aggregation / polymerization of units of the soluble pool to grow the insoluble pool and the disassembly / solubilization / depolymerization of the insoluble pool into units of soluble pool. It is assumed that the assembly process is a function of the soluble pool only, whereas the disassembly process is a function of the insoluble pool only. The simplest case is to consider a linear model that assumes linear exchanges between the two pools. The linear model can be stated as follows: where k (⋅) (resp. k (⋅)) is the rate of assembly of the soluble pool (resp. disassembly of the insoluble pool). Both rates can either be constant, k and k , or space-dependent k (x) and k (x).

In the second case, the turnover is assumed to depend on the enzymatic activities [22]. For instance, the solubilization of the insoluble pool into soluble proteins (disassembly) is triggered by a kinase activity and the assembly of insoluble pool is induced by the dephosphorylation of soluble proteins by a phosphatase [23]. The turnover term is assumed to be of the Michaelis-Menten form stated as: where k (⋅) (resp. k (⋅)) is the maximum rate of assembly of the soluble pool (resp. disassembly of insoluble pool) and k (resp. k ) is the concentration at which the assembly (resp. disassembly) rate is half of k (⋅) (resp. k (⋅)). When Michealis-Menten dynamics are used, the model is called nonlinear. Both rates can either be constant or space-dependent describing the intracellular localization of the post-translational modification enzymes.

Profiles of assembly and disassembly rates

As previously mentioned, the assembly and disassembly rates, k (⋅) and k (⋅), used in the linear and nonlinear models can be constant or space-dependent functions. The profile (shape) of the space-dependent function k (x) is derived from the spatial profile of regions of assembly (Sources) measured in [12] (Fig 2.3). Details of the derivation of the Sources k (x) are given in Appendix 2.

Two types of shapes for k (x) are assumed to represent two types of localization of the disassembly in cells. First, similarly to the assembly rate, the profile of the disassembly rate is deduced from the experimental data published in [12] (Fig 2.3); the spatial profile of the disassembly regions, Sinks, is used to build the shape of the first space-dependent disassembly rate. This disassembly rate is called of type Sinks. Second, it is assumed that disassembly is localized around the nucleus; a mollified step-function is designed to describe this assumption. This second space-dependent disassembly rate is called of type Mollify. Details on the derivation of the two k (x) rates are given in Appendix 3.

For the sake of illustration, the two profiles of k (⋅) and the three profiles of k (⋅) are shown in Fig 4.

Fig 4

Possible profiles for the assembly and disassembly rates.

4.1: k = 10−3 s −1 when the linear model is used (resp. μM.s −1 when the nonlinear model is used) and k (x) is computed using (14). 4.2: k = 10−3 s −1 when the linear model is used (resp. μM.s −1 when the nonlinear model is used) and k (x) are computed using (15) for Sinks and (16) for Mollify, respectively. Details on the derivation of k (x) and k (x) rates are given in Appendix 2 and Appendix 3, respectively. Parameter values in (14), (15) and (16) are chosen in such a way that their profiles give the same total amount of assembly / disassembly over the spatial domain Ω.

Accounting for the three modes of transport and considering the turnover described by either linear and nonlinear models with 6 possible combinations of the profiles for the assembly and disassembly rates, a collection of 36 scenarios (models) is defined; each scenario follows the general form stated by system (1). Details of the 36 scenarios are given in Fig 5. All scenarios are considered with the same initial conditions given by where f 0(x), defined in (17), is chosen to be a mollified version of the profile of the assembled keratin at time t 0 = 24 hours averaged over all the normalized cells (Fig 6). Details of the derivation of f 0(x) are given in Appendix 4. Initial conditions describe the observed fact that the soluble pool (resp. insoluble pool) represents 5% (resp. 95%) of the total keratin material [24]. All scenarios are considered with boundary conditions describing the impermeability of the plasma membrane for the keratin material where J is the flux of the soluble pool and J is the flux of the insoluble pool as defined below system (1). The model parameters used in all the scenarios are listed in Table 1.

Fig 5

The 36 scenarios to be considered.

Top: Number in parentheses is the scenario index i. Bottom: The numerical value 1 denotes that the process of interest is in Scenario i.

Fig 6

Initial profile of the insoluble pool f 0(x) defined in (17).

Comparison between mathematical models and experimental data

Parameter estimation

Let p denote the set of all the model parameters for each scenario. To estimate the optimal set of parameter values p for each scenario, the solution of each scenario is compared to experimental data using the following objective function (an estimation of the distance between the experimental data and the model response): where t is equal to 48 hours. The experimental data at t is represented by f (x) and the model solution at t of the considered scenario evaluated with the parameter values p is I(x, t , p). To obtain the model solution for each scenario, the corresponding system is numerically integrated from t 0 = 24 hours to 50 hours using the MATLAB solver for partial differential equations, pdepe [25]. Solutions are computed at N locations of the spatial domain Ω with N = 200. To carry out computations, the raw data in which the concentration of the assembled keratin material is known at 623 spatial locations is approximated by f (x) defined in Appendix 4 by (18); f (x) is the fit of the average profile of the assembled keratin measured after 48 hours on the normalized cells (see the black profile in Fig 2.1).

The estimation of the parameter values for the 36 scenarios, i.e. the determination of parameter values that provide the best fit to the experimental data, is done by minimizing the objective function, Ψ(p), such that Minimization of the objective function is done by a parallelizable genetic algorithm described in [26].

In this study, only constant parameter values of the turnover reactions are optimized. When the linear model is considered, for each scenario, two parameters k and k are estimated. When the nonlinear model is considered, four parameters, k , k , k and k , are estimated for each scenario. In all scenarios, the diffusion coefficients for the soluble and insoluble pools as well as the active transport speeds v(x) and u(x) are considered as fixed/known parameters or functions and are set to values measured in previous studies [9, 11, 12] (Table 1). The diffusion coefficients are taken as D = 0.88μm 2/s and D = 9.5 × 10−4 D , and the “constant speed” u in v(x) as defined in (12) is set to be u = 0.0025μm/s [9, 18].

Model selection

In order to select the best model out of the 36 scenarios, we use an information-theoretic approach. Specifically, we employ the Akaike information criterion (AIC) [27, 28] to select from all the scenarios, the best model that best captures experimental observations given the collection of models considered in this study. Since we use the Least Squares principle to estimate the values of the constant free parameters (2 free parameters for linear models and 4 for nonlinear models), the Akaike criterion for Scenario i, AIC , takes the following form: where is the estimate of the variance with the residual (9) estimated for Scenario i and N the size of the sample (N = 200). K is the number of estimated parameters; i.e. the number of free parameters in Scenario i plus one for the estimate of the variance (K = 3 for linear model and K = 5 for nonlinear models). The scenario with the lowest AIC value is the best model. The AIC selects a model with the least number of parameters that best fits experimental observations. To rank and compare scenarios the Akaike weights w are calculated and these are known as the weights of evidence in favor of Scenario i being the actual best model given the experimental data and the collection of scenarios considered. The Akaike weights are expressed as follows: where R represents the number of scenarios considered (R = 36) and Δ being the difference in AIC with respect to the AIC of the preferred scenario min AIC = AIC . It must be noted that the Akaike weights w sum to 1 and are interpreted as the probability that Scenario i is the best model given the experimental data and the collection of scenarios considered. The models, ranked from the largest to the smallest Akaike weights, whose Akaike weights sum to 0.95 form the confidence set of the models that captures, more faithfully, experimental data. The ratio of the Akaike weights w /w (also known as the evidence ratio) is used to compare model pairs. Furthermore, the relative importance of a process can be estimated by summing the Akaike weights of all scenarios involving the process of interest. We will denote by the sum of the Akaike weights of the scenarios including the process of type *. This sum can be interpreted as the probability that the process of type * is the best type of the process given the experimental data and the collection of models considered.

Results

Best scenario

We employ a two-step process in order to find the best scenario. First, the best fit for each scenario is found by minimizing the objective function (9). Secondly, the Akaike information criterion (10) is used to select the best of the best fits out of all the scenarios.

The best fit for each of the 36 scenarios is presented in Fig 7 in the following order:

Fig 7

Best fit for each of the 36 scenarios.

Red curve is the response at 48 hours of Scenario i. Black circles are data at 48 hours as in Fig 2.1.

Row-wise:

Scenarios in the first three rows (1 to 3) have linear terms for turnover.

Scenarios in the last three rows (4 to 6) have Michaelis-Menten type turnover terms.

Rows 1 and 4 are scenarios with constant disassembly rate.

Rows 2 and 5 are scenarios with localized disassembly rate of type Sinks.

Rows 3 and 6 have scenarios with localized disassembly rate of type Mollify.

Column-wise:

The first two columns (1 and 2) include scenarios with diffusion alone without drift.

Columns 3 and 4 display scenarios with drift with an almost constant speed.

The last 2 columns (5 and 6) are scenarios with drift with variable speed.

Odd columns (1, 3 and 5) are scenarios with constant assembly rate.

Even columns (2, 4 and 6) are scenarios with localized assembly rate of type Sources.

According to the AIC values, the best model, the best of the best fits, is identified as being Scenario 21, AIC 21 = min AIC (3 column of Table 2). Since Scenario 21’s Akaike weight is over 0.95 (5 column of Table 2), it is the only model out of the set considered that satisfies the confidence criterion; Scenario 21 matches very well experimental observations and data. The second best model is Scenario 31. The evidence ratio of scenarios 21 and 31 w 21/w 31 is equal to 36.67 × 106; i.e. Scenario 21 is about 36 millions times more likely than Scenario 31 to be the best model given the experimental data and the collection of models considered. We consider this to be strong evidence in support of Scenario 21.

Table 2

Results of the model selection for the best fit of each of the 36 scenarios.

Scenario i	K	AICi	Δi	wi	Rank
1	3	2012.996	510.3519	0	13
2	3	2019.144	516.4999	0	19
3	3	2018.634	515.9895	0	18
4	3	2512.249	1009.604	0	36
5	3	2283.252	780.6081	0	28
6	3	2340.069	837.4254	0	31
7	5	2011.436	508.7923	0	11
8	5	2011.724	509.0804	0	12
9	5	2013.609	510.9647	0	14
10	5	2014.842	512.1975	0	15
11	5	2015.148	512.5043	0	16
12	5	2017.034	514.3897	0	17
13	3	1839.930	337.2861	0	7
14	3	2479.026	976.3817	0	35
15	3	1738.738	236.0939	0	5
16	3	2114.919	612.2750	0	23
17	3	2458.497	955.8532	0	34
18	3	2002.972	500.3279	0	9
19	5	1668.771	166.1274	0	3
20	5	2387.721	885.0770	0	32
21	5	1502.644	0	0.99999	1
22	5	2073.809	571.1653	0	21
23	5	2338.356	835.7119	0	29
24	5	1893.875	391.2309	0	8
25	3	1703.773	201.1292	0	4
26	3	2410.565	907.9213	0	33
27	3	2006.784	504.1403	0	10
28	3	2067.488	564.8442	0	20
29	3	2243.066	740.4220	0	27
30	3	2145.639	642.9953	0	26
31	5	1537.479	34.83477	2.7 × 10−8	2
32	5	2339.200	836.5561	0	30
33	5	2119.447	616.8033	0	24
34	5	2074.025	571.3811	0	22
35	5	1834.186	331.5420	0	6
36	5	2143.595	640.9505	0	25

The numerical value 0 in the 5 column denotes a value lower than 10−12. In the 6 column, “Rank” denotes the ranking of scenarios, i.e. the descending order of Akaike weights w . The use of weights allows the “quantitative” comparison of the adequacy of scenarios because of the normalization to 1.

Scenario 21 is characterized by an inward drift with an almost constant speed for the insoluble pool, turnover terms of Michaelis-Menten type, a constant assembly rate and a disassembly rate of type Mollify. The profiles of assembly and disassembly rates are displayed in Fig 8. For this scenario, the estimated optimal parameter values are k = 9.3819μM/s, k = 570.73μM, k = 0.9998μM/s leading to k = 1.976μM/s in (16) and k = 976.07μM. Based on the Michaelis-Menten constants for the assembly and disassembly processes, since k < k , an enzyme that would be involved in the solubilization of the assembled keratin material requires a higher substrate concentration to achieve a given reaction speed than an enzyme that would be involved in the assembly of soluble proteins. Snapshots of the soluble and insoluble pool profiles taken every 30 minutes from t 0 to t are displayed in Fig 8. After only 2 hours, the solution stabilizes to its final profile. Scenario 21 preserves the repartition of the keratin material between the soluble and insoluble pools over time, about 95% of the keratin material is assembled to form the insoluble pool. Finally, the characteristic time scales of the passive transport (τ ), active transport (τ ) and turnover (τ ) for Scenario 21 are estimated by using an adimensionalization of system (1) corresponding to Scenario 21. Details on calculations are given in Appendix 5. The time scales of the processes included in Scenario 21 are ordered as follows: Using Peclet’s number (Pe = τ /τ ≫ 1), it is found that the transport of the assembled keratin material by drift is faster than by diffusion, the dominant mode of transport is the active transport. Using Damköhler number (Da = τ /τ ≫ 1), it is found that the active transport time scale is greater than the reactive time scale; overall, Scenario 21 is controlled by active transport.

Fig 8

Details about the best model Scenario 21.

8.1: Snapshots from 24 hours (gray) to 48 hours (black) taken every 30 minutes of the profile of the soluble pool over the cell. 8.2: Snapshots from 24 hours (gray) to 48 hours (black) taken every 30 minutes of the profile of the insoluble pool over the cell. 8.3: The distribution of the keratin material between the soluble and insoluble pools over time. 8.4: The profiles of assembly and disassembly maximal rates used with the Michaelis-Menten type turnover terms in Scenario 21. The estimates are k = 9.3819μM/s, k = 570.73μM, k = 0.9998μM/s leading to k = 1.976μM/s in (16) and k = 976.07μM/s.

The complete ranking based on the Akaike weights ω of all the scenarios is given in the 6 column of Table 2. It is worthwhile remarking that differences in AIC , Δ, could have been used for ranking purposes; in this case, we would have obtained the same ranking. The rationale of using ω is that it allows us to quantify how preferably each candidate model is via the normalization to 1.

Importance of the type of process

Using the sums of the Akaike weights, , we investigate several questions to evaluate the relative importance of the different types considered for a process (Tables 3 and 4). The set of the 36 scenarios is partitioned into categories with respect to the types of processes considered. For instance, considering the transport process, the collection of the 36 scenarios is partitioned into 3 categories (Fig 5): scenarios with no drift (Scenarios from 1 to 12), scenarios with an almost constant speed drift (Scenarios from 13 to 24) and scenarios with a variable speed drift (Scenarios from 25 to 36). The sum of the Akaike weights of each category is then calculated, compared and ordered to characterize which type is more likely to be present. In what follows the sign > denotes the relative importance, measured in probabilistic terms, of the type of process.

Table 3

Importance of the type of process.

Type of Transport			Type of Turnover Term		Type of Assembly		Type of Disassembly
wNoDrift+=0	wCstDrift+=1	wVarDrift+=2.7×10^−8	wLinear+=0	wNonlinear+=1	wCstAss+=1	wSource+=0	wCstDis+=2.72×10^−8	wSink+=0	wMoll+=1

denotes the sum of the Akaike weights of scenarios including the type * of the process. Type of transport—No drift for the insoluble pool (No drift) vs Drift with almost constant speed (Cst drift) vs Drift with variable speed (Var. drift). Type of turnover terms—Linear vs Nonlinear. Type of assembly—Non-localized (constant) vs Localized of type Sources. Type of disassembly—Non-localized (constant) vs Localized of type Sinks vs Localized of type Mollify.

Table 4

Importance of the type of process.

Scenario i	Cst Ass. Dis.	Cst Ass. Sinks	Cst Ass. Moll.	Sources Cst Dis.	Sources Sinks	Sources Moll.
1	1
2		1
3			1
4				1
5					1
6						1
7	1
8		1
9			1
10				1
11					1
12						1
13	1
14		1
15			1
16				1
17					1
18						1
19	1
20		1
21			1
22				1
23					1
24						1
25	1
26		1
27			1
28				1
29					1
30						1
31	1
32		1
33			1
34				1
35					1
36						1
	wCstAssDis+=2.7×10^−8	wCstAssSink+=0	wCstAssMoll+=1	wSourceDis+=0	wSourceSink+=0	wSourceMoll+=0

Combinations of assembly/disassembly rate profiles—Constant assembly and disassembly rates vs constant assembly rate and disassembly rate of type Sinks vs constant assembly rate and disassembly rate of type Mollify vs assembly rate of type Sources and constant disassembly rate vs assembly rate of type Sources and disassembly rate of type Sinks vs assembly rate of type Sources and disassembly rate of type Mollify. Top: The numerical value 1 denotes that the combination of processes of interest is in Scenario i. Bottom: denotes the sum of Akaike weights of scenarios including the process combination as indicated at the top.

What is the most likely type of transport for the keratin material in cells given the experimental data and the model collection considered? Is it that there is no drift of the insoluble pool (diffusion only for both pools), drift with an almost constant speed for the insoluble pool or drift with a variable speed for the insoluble pool? Each category includes 12 scenarios (see 3 to 5 columns in Fig 5). From the Akaike weights, it is obvious that scenarios including diffusion only (with no drift) are not supported by experimental data since (Table 3). Diffusion alone is not enough to explain the experimental data. Given the data and the model collection, the type of transport can be ordered as follows (Table 3):

What is the most likely type of turnover between the soluble and insoluble pools given the experimental data and the model collection considered? Is linear exchange or Michaelis-Menten type (underlying an enzyme activity) more likely? Each category includes 18 scenarios (see 6 to 7 columns in Fig 5). From Table 3, enzymatic activities are more likely to be present. Hence

What is the most likely rate profile of the assembly process of the keratin material in cells given the experimental data and the model collection considered? Is a constant assembly rate or assembly mainly localized at the cell membrane periphery more likely? Each category includes 18 scenarios (see 8 to 9 columns in Fig 5). From Table 3, the non-localized rate of assembly is preferred to the rate of the localized assembly; . Hence

What is the most likely rate profile of the disassembly process of the keratin material in cells given the experimental data and the model collection considered? Is a constant disassembly, disassembly of type Sinks or disassembly localized around the nucleus more likely? Each category is composed of 12 scenarios (see 10 to 12 columns in Fig 5). From Table 3, none of the scenarios that include the disassembly rate of type Sinks is supported by experimental data.

What is the most likely combination of the assembly and disassembly rate profiles given the experimental data and the model collection considered? Is a combination of constant assembly and disassembly rates, a constant assembly rate and disassembly rate of type Sinks, a constant assembly rate and disassembly rate of type Mollify, an assembly rate of type Sources and constant disassembly rate, an assembly rate of type Sources and disassembly rate of type Sinks or an assembly rate of type Sources and disassembly rate of type Mollify more likely? Now, interactions between the assembly and disassembly processes are investigated. Six categories are defined. In each category, there are 6 scenarios as listed in the top of Table 4. Overall, as in Scenario 21, the combination of non-localized assembly and disassembly localized around the nucleus is preferred (Table 4).

Discussion

The primary aim of the present work is to investigate, through mathematical modeling driven by experimental observations, mechanisms contributing to the organization of the keratin material in epithelial cells and subsequently to test the biological model for the keratin dynamics proposed by Leube and Windoffer et al. in 2011 [10, 11]. From first principles, we formulate a collection of mathematical models capturing various combinations of biological processes describing the spatio-temporal dynamics of the keratin material in epithelial cells. By using techniques in parameter estimation, we find optimal reaction kinetic parameter values for each model that best captures experimental observations. We go one step further and employ an information-theoretic approach for model selection to determine the best model among all the best fit models that captures the key processes of the experimental data.

As previously highlighted, the model framework and the description of processes considered are driven by experimental data and observations. Experimental data used in this study provide the spatial distribution of the assembled keratin concentration at 24 and 48 hours. As only macroscopic information on the keratin organization is available, the models only describe the keratin material as soluble or assembled and consider exchanges between these soluble and insoluble pools combined with transport events. For instance, as experimental data identify the existence of regions with preferential assembly and disassembly (Fig 2.3), the assembly and/or disassembly processes are assumed to be localized in some scenarios. When observing the perpetual inward motion of the keratin network (see S1 Video), the “natural assumption” for the transport of the assembled keratin would be the existence of an inward active transport. However, we decide to consider scenarios with only passive transport (no active transport). Why? First, it is well known that combining appropriate reaction terms and diffusion (only) can lead to the emergence of complex behavior such as traveling waves and/or pattern formation (see some examples in [29]). Secondly, to reinforce the predictive power of our study. In our approach, we do not only design mathematical models but also calibrate models (parameter estimation) and then evaluate (model selection) how each model performs. To highlight and confirm the importance or existence of a given process, scenarios in which the process is absent must also be considered and evaluated. Our strategy of considering “unrealistic scenarios” and “realistic scenarios” in the collection of models to be evaluated mimics the biological experimental protocols such as knockout and controls. For instance, it has been shown that the keratin assembly / disassembly depend on post-translational modifications of keratins due to enzymatic activities [23]. In our study, a total of 36 scenarios that divide into 18 scenarios with a turnover of type linear (no enzymatic activities) and 18 scenarios with a turnover of type nonlinear (enzymatic activities) are investigated. There is a one-to-one correspondence between the 18 scenarios of the 2 groups. When model performances (Akaike weights) are compared, scenarios having a nonlinear turnover perform better than the corresponding ones with a linear turnover. Hence, the enzymatic activity is detected by the model selection as existing but also preferable. This outcome allows us to judge that our approach combining the mathematical modeling and model selection performs correctly, it gives us the confidence and trust in our conclusions.

Following this methodology, it turns out that Scenario 21 was the best of the best fits. Scenario 21 is in good agreement with the biological model proposed in [10, 11]. Scenario 21 and the biological model share the following common key features: diffusion of the soluble pool, inward active transport of the insoluble pool and disassembly of the insoluble pool localized around the nucleus. Moreover, as experimentally observed, Scenario 21 preserves the repartition of the keratin material over time; the keratin material is mainly insoluble in epithelial cells. Both the proposed biological model and Scenario 21 hypothesize an inward active transport for the assembled keratin material. Examining the Akaike weights of the models with no active transport (Scenarios 1 to 12 in Table 2, 1 and 2 columns in Fig 7 and 1 to 3 columns in Table 3), we go one step further and show that the active transport is a requirement to explain the experimental data and that the experimentally observed inward motion of the assembled keratin is not an emergent phenomenon but it is due to an active transport. Furthermore, comparing the characteristic time scales of processes, we establish that the keratin dynamics are mainly controlled by the active transport in Scenario 21. Interestingly, Scenario 21 concurs with the proposed biological model about the existence of a disassembly localized in the perinuclear region. An experimental protocol must now be developed to further work out the details of this conclusion. It is worth pointing out that some characteristics of the proposed biological model were not tested due to the form of the mathematical models. For instance, the scenarios considered here describe the turnover in terms of the assembly of soluble proteins and disassembly of aggregated proteins. In the biological model, the nucleation of filaments, i.e. the initiation of filaments, is explicitly described and localized at the cell periphery. However, the nucleation process is not explicitly described in any of the present models. It is worth mentioning that at the beginning of this work a nucleation term was included in the models; the effect of this term did not change significantly the results; hence, to reduce the complexity of models and the number of parameters the nucleation term was subsequently dropped. Investigation, by mathematical modeling, of the nucleation process is currently being carried out in a separate and on-going study.

When only scenarios with a variable drift are considered (Scenarios 25 to 36, 5 and 6 columns in Fig 7), the preferred combination of assembly and disassembly rate profiles is a constant rate of assembly and disassembly as in Scenario 31, which is the second best model. The best model, Scenario 21, and the second best model, Scenario 31, share the non-compartmentalization of the assembly process. However, in Scenario 21 the active transport has an almost constant speed whereas in Scenario 31 the speed is variable. For the disassembly, this process is non-localized in Scenario 31 whereas it is localized at the perinuclear region in Scenario 21. The spatial variability of the speed u(x) in Scenario 31, in particular, the almost-zero speed at the nucleus locations (see for x ∈ [−7.5, 7.5] in Fig 3) “compensates” for the non-compartmentalization of the disassembly process. After comparing the features of Scenarios 21 and 31, we inspect their profiles. The profile obtained with Scenario 21 fits very well the experimental data on non-perinuclear locations (see for x ∉ [−7.5, 7.5] in Fig 7.33). On the other hand, the profile resulting from Scenario 31 fits very well experimental data on perinuclear locations (see for x ∈ [−7.5, 7.5] in Fig 7.23). As an improvement of Scenario 21, we would expect that a wider decay in speeds around the nucleus modeled, for instance, with a smaller value of a in (12) would result in a better fit of the perinuclear region.

Scenarios 29 and 35 include the variable speed measured from experimental data, a localized assembly rate having the spatial profile deduced from Sources and a localized disassembly rate whose shape is derived from the profile of the Sinks. All the characteristics extracted from the experimental data are included in Scenario 29 (linear model) and Scenario 35 (nonlinear model). If one wanted the best model that takes into account all the experimental features, then Scenario 29 or 35 would be the best scenario. Scenario 29 is ranked 27 and Scenario 35 is ranked 6 (Table 2). Similarly to the general trend, for the same set of assumptions, the use of the Michaelis-Menten type turnover term generally provides a better representation of experimental data than the use of the linear turnover term. The evidence ratio of Scenarios 35 and 29 ω 35/ω 29 is ridiculously large; Scenario 35 is much more adequate to represent the experimental data than Scenario 29. However, Scenario 35 is still only ranked 6. The failure/mismatch of Scenario 35 might be explained by the redundancy of the information existing in the variable speed and the net assembly/disassembly region profiles extracted from the experimental data. Furthermore, according to the Akaike weights, disassembly rates of type Sinks (the type deduced from experimental data) are less likely to occur given the experimental data and the collection of models considered. A similar conclusion is obtained for the assembly rate of type Sources. This could be a consequence of our too conservative interpretation of the regions of preferential assembly or disassembly. In our assumptions, in zones of preferential assembly, the disassembly rate is set to be about null and vice versa (Fig 2.3).

In summary, to model the keratin dynamics in epithelial cells we characterized the keratin material into two pools, the soluble and the insoluble pool; and, the events considered are the turnover and transport for both pools. The modeling assumptions used, for instance, the diffusion of the soluble pool, are based on biological observations and experimental data. The collection of the models considered in this study is designed to answer a set of questions such as “what is the mode of transport of the keratin material in epithelial cells?”. After optimizing parameter values, model selection and evaluation methods applicable to non-nested models are used to discriminate between the candidate models, identify the best model and quantify how models under consideration are adequate to explain the experimental data. Note that the ranking of the models (scenarios) and the relative importance of the different types of processes (Tables 2–4) are only valid in the context of the experimental data and the set of the candidate models considered here. For instance, considering other hypotheses such as the non-negligence of the anterograde motion of the assembled keratin material or the stabilization (protection against disassembly) of the keratin filaments involved in the nuclear cage would have led to a different collection of scenarios in which our best scenario could have failed to be the best one. Furthermore, as some modeling assumptions are directly derived from experimental data, missing information in experimental data might have been prejudicial for the correct approximation / modeling of processes; for instance, missing information about the speed of assembled keratin in the cell periphery result in small values for the variable speed close to the plasma membrane (Fig 2.2). Keeping in mind the limitations of our approach, important conclusions have been reached such as

an active transport of the assembled keratin material is required, thereby confirming recent experimental observations [10, 11],

enzymatic activities regulating the assembly / disassembly are more likely to occur,

the assembly process is more likely to be non-compartmentalized in cells,

last but not least, a unique best model strongly supported by experimental data is identified, this scenario is in good agreement with the biological model previously proposed in [10, 11]. Interestingly, the best scenario supports the perinuclear localization of the disassembly process hypothesized in the biological model.

Mathematical models of the keratin intermediate filament organization in cells were previously proposed [30-34]; however, none of these models included the effects of transport of the assembled material on its organization. More importantly, only the behavior of models were validated qualitatively, no comparisons to experimental data were carried out. It is worth noting that other studies of the intermediate filaments dynamics combining mathematical modeling and experimental data approaches were carried out but on neurofilaments in neurons (see for example [35, 36]).

Conclusion

Given the experimental data published in [12], through modeling and simulations, we investigate the effects of the interplay between turnover and transport on the keratin spatio-temporal organization in epithelial cells. Out of all the scenarios investigated, a scenario strongly supported by experimental data is found that best captures most of the hallmarks of the experimental observations. This scenario predicts the diffusion of soluble keratin, an inward active transport of the assembled keratin and the disassembly localized around the nucleus triggered by enzymatic activities as well as the assembly process that is non-compartmentalized over the cell. The value of our models is reflected in their predictive nature; first, the approach predicts the localized disassembly at the perinuclear region and second, that the experimentally observed inward motion is not an emergent behavior but that it is an inherent property of the organization of the keratin material in epithelial cells and it is due to an active transport.

Appendix 1

Based on experimental observations, the speed of the assembled keratin material is assumed to decay to almost zero around the nucleus. The nuclear envelope positions are at x = ±7.5μm. In order to describe the decay of the speed around the nucleus such that it is almost constant, the following function is used: with a = 0.05 and u being in the range given in Table 1. For numerical simulations, u is set to 0.0025μm/s which represents the average value of speeds measured in [12].

For the variable speed case, a symmetrical function over the spatial domain Ω is sought to describe the speed u(x). Averaging over the symmetrical spatial locations values of the average speed measured in [12] (Fig 2.2) and curve fitting with a sum of Gaussian functions of these values, an estimate of u(x) is obtained and expressed as follows: with the following coefficients: a 1 = 0.003372, b 1 = 17.39, c 1 = 7.577, a 2 = 0.003378, b 2 = −17.41 and c 2 = 7.546. The estimate u(x) is almost zero at the nucleus-locations and is symmetrical around the center of the cell (Fig 9).

Fig 9

Estimate of the space-dependent speed.

Estimate of the space-dependent speed given in (13) for the variable speed u(x) of the active transport of the insoluble pool.

For numerical simulations, when needed, the following function with a = 1 is used as a “smooth analogue” of sgn(x).

Appendix 2

To obtain the space-dependent function k (x) the profile of the regions of assembly (Sources) published in [12] (Fig 2.3) is first made symmetrical by averaging values at the symmetrical spatial locations and then the symmetrical Sources profile is fitted using the sum of Gaussian functions: The best fit is obtained with the following coefficients: a 1 = 1.684 × 105, b 1 = 20.13, c 1 = 2.105, a 2 = 8.493 × 104, b 2 = −19.04, c 2 = 1.374, a 3 = 1.333 × 105, b 3 = −20.74, and c 3 = 1.734. The value of the coefficient k determines the maximal value of the peaks. When with k = k × 10−2, the total amount of assembly over the cell is the same as that of the case of a constant assembly of level k . The value k in k is determined by fitting model solutions to experimental data. An illustration of the shape of k (x) is given in Fig 10.

Fig 10

Profile for localized assembly rate of type Sources.

The profile for k (x) defined in (14) with k = 1, obtained by fitting the profile of assembly regions Sources published in [12] and shown in Fig 2.3.

Appendix 3

To obtain the first space-dependent function k (x), the profile of regions of disassembly (Sinks) published in [12] (Fig 2.3) is made symmetrical and fitted using the sum of Gaussian functions: The best fit of the symmetrical profile of Sinks is obtained with the following coefficients: a 1 = −3.383 × 105, b 1 = −11.51, c 1 = 4.158, a 2 = −1.004 × 105, b 2 = −14.03, c 2 = 2.236, a 3 = 2.445 × 105, b 3 = −11.39, c 3 = 4.339, a 4 = −2.032 × 104, b 4 = 5.922, c 4 = 5.385, a 5 = −2.058 × 104, b 5 = −6.077, c 5 = 5.778, a 6 = 3.793 × 105, b 6 = 11.37, c 6 = 4.322, a 7 = −1.008 × 105, b 7 = 14.03, c 7 = 2.238, a 8 = −4.742 × 105, b 8 = 11.46, and c 8 = 4.196. The coefficient k determines the amplitude of the peaks. When , the total amount of disassembly over the cell is the same as that of the case of the constant function of level k . The value k in k is determined by fitting model solutions to experimental data. An illustration of the shape of k (x) obtained from Sinks is given in Fig 11.

Fig 11

Profile for localized disassembly rate of type Sinks.

The profile for k (x) defined in (15) and k = 1, obtained by fitting the profile of disassembly regions Sinks published in [12] and shown in Fig 2.3.

A second function is hypothesized to represent the localized disassembly k (x) around the nucleus. A mollified piecewise function is used and is of the form where a 1 = −15μm, a 2 = −7.5μm, a 3 = 7.5μm, a 4 = 15μm and ϵ = 1. When k = 2(k −k ) with k = k × 10−2, the total amount of disassembly over the cell is the same as that of the case of a constant disassembly of level k ; k is determined by fitting the model solutions to experimental data. An illustration of the shape of (16) is given in Fig 4.

Appendix 4

To define the initial condition f 0(x), the experimental profile [12] of the assembled keratin material measured at 24 hours is used (gray circles in Fig 2.1). The polynomial for which p 8 = 5.441 × 10−8, p 7 = 3.397 × 10−21, p 6 = −5.379 × 10−5, p 5 = −2.077 × 10−18, p 4 = 0.01062, p 3 = 2.801 × 10−16 p 2 = 0.4104, p 1 = 1.4 × 10−14 and p 0 = 506.5 fits well the data at 24 hours (gray curve in Fig 2.1). The function P(x) is then modified to obtain a function satisfying the boundary conditions defined in (7). The initial condition f 0(x) describing the profile of the insoluble pool at 24 hours is expressed as follows: with a 1 = −21μm, a 2 = −a 1, ϵ = 1, c = 50μM and . The function f 0(x) is graphed in Fig 6.

The average profile of the assembled keratin material computed on the normalized cells after 48 hours of seeding [12] is represented by: where p 4 = −0.003255, p 3 = 2.61 × 10−17, p 2 = 0.4899, p 1 = 1.558 × 10−15 and p 0 = 604.1. f (x) is graphed in Fig 2.1 (black curve).

Appendix 5

The system corresponding to Scenario 21 is now nondimensionalized to estimate the characteristic time scale of each of the processes involved. Consider new independent variables z and τ and new dependent variables and defined such that: Scenario 21 takes then the following form: with f(Az) being the functional component of v(x) (12) defined as f(Az) = (1 − exp(−a(Az)2)) with a = 0.05 and k (Az) being the function k (⋅) defined in (16) and evaluated at Az. Taking A = ℓ (ℓ = 2L), B = ℓ/u, s = k and i = k then and the adimensional version of the system is The constant is the Péclet number of the insoluble pool that compares active transport and diffusion processes. The characteristic length ℓ is chosen to be equal to the length of the spatial domain Ω, ℓ = 2L. The Damköhler number compares reaction and active transport processes. The characteristic active transport time is τ = ℓ/u as defined in the Péclet number. The characteristic reaction time is τ = 1/κ for which . Note that sgn(z)f(Az) is about ±1.

Recall the parameter values obtained from the parameter estimation: k = 9.3819μM/s, k = 570.73μM, k = 976.07μM/s and k = 1.976μM/s in (16) that correspond to a constant disassembly rate of level k = 0.9998μM/s. Hence, from the definition of k , k (Az) can be approximated by k = 0.9998. The values of the fixed parameters are D = 0.88μm 2/s, D = ϵD with ϵ = 9.5 × 10−4 and u = 0.0025μm/s. The following estimates are then obtained:

Adimensional parameters: k (Az)/k = k /k = 0.1 and k /k = 1.7;

Diffusion time scale: τ = ℓ2/D.

– For the soluble pool: ,

– For the insoluble pool: ;

Active transport time scale: τ = ℓ/u = 1.8 × 104 s (time the fluid needs to flow through the characteristic length);

Reaction time scale: k /k ≈ 61s (time for reaction to equilibrate).

To determine what mode of transport (passive vs active) and what process (transport vs reaction) dominate the dynamics of the keratin material, the Péclet and Damköhler numbers are used:

Pe = (uℓ)/(ϵD ) ≈ 135 ≫ 1. Diffusion is small compared to active transport for the insoluble pool, the transport of the assembled keratin material by drift is faster than by diffusion, active transport is the dominant mode of transport;

Da = (ℓk )/(uk ) ≈ 296 ≫ 1. Active transport time scale is greater than the reaction time scale; the overall process in controlled by active transport.

Dynamics of the keratin network in a cell.

Time-lapse fluorescence microscopy of hepatocellular carcinoma-derived PLC clone PK18-5 stably expressing fluorescent fusion protein HK18-YFP [8] depicting the dynamic properties of the keratin filaments over a time period of 15 hours. Bar 10 μm.

(AVI)

Click here for additional data file.