Keratin Dynamics: Modeling the Interplay between Turnover and Transport
Stephanie Portet, Anotida Madzvamuse, Andy Chung, Rudolf E. Leube, Reinhard Windoffer
KKeratin Dynamics: Modeling the Interplay betweenTurnover and Transport
St´ephanie Portet , Anotida Madzvamuse , Andy Chung , Rudolf E. Leube , ReinhardWindoffer Abstract
Keratin are among the most abundant proteins in epithelial cells. Functions of thekeratin network in cells are shaped by their dynamical organization. Using a collectionof experimentally-driven mathematical models, different hypotheses for the turnoverand transport of the keratin material in epithelial cells are tested. The interplaybetween turnover and transport and their effects on the keratin organization in cells arehence investigated by combining mathematical modeling and experimental data.Amongst the collection of mathematical models considered, a best model stronglysupported by experimental data is identified. Fundamental to this approach is the factthat optimal parameter values associated with the best fit for each model areestablished. The best candidate among the best fits is characterized by the disassemblyof the assembled keratin material in the perinuclear region and an active transport ofthe assembled keratin. Our study shows that an active transport of the assembledkeratin 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 mechanicalstress resistance of epithelial cells and contributes significantly to epithelialstiffness [26, 28]. The importance of keratins for epithelial tissue stability is reflected bya large group of genetic skin blistering diseases that are caused by point mutations inkeratin-encoding genes [6, 11]. The mechanical functions not only rely on staticresilience but also necessitate a high degree of plasticity, for example in migrating cellsduring wound healing [5]. The current view is that keratins act as general stressabsorbers protecting epithelial cells not only against mechanical insults but also againstirradiation or osmotic and microbial challenges. Thus, keratins are involved in heatshock response, apoptosis and organelle homeostasis [32]. Furthermore, functionsaffecting processes such as proliferation, differentiation and inflammation are alsodependent on keratins (see recent reviews in [11, 23]). 1/29 a r X i v : . [ q - b i o . S C ] A p r igure 1. Keratin network. Image taken from a time-lapse fluorescence recording ofa single hepatocellular carcinoma-derived PLC cell of clone PK18-5 [30] stablyexpressing fluorescent fusion protein HK18-YFP consisting of human keratin 18 andenhanced yellow fluorescent protein. Bar 10 µm .All of these functions are tightly coupled to the keratin network dynamics (see S1Video). Examination of cultured epithelial cells producing fluorescent keratins hasprovided evidence for the different mechanisms that are involved in the continuousrenewal and reshaping of the keratin system [16]. On the basis of these observations, abiological model of the keratin cycle has been proposed by Leube and Windoffer et al. in [17, 33]. This biological model takes into account the assembly/disassembly andtransport of keratins. In this study, it was proposed that assembly and disassemblyoccur in topologically defined regions with assembly taking place predominantly in thecell periphery while disassembly takes place primarily in the perinuclear region. Thebiological model further postulates active transport of insoluble assembly stages ofkeratins toward the nucleus and rapid diffusion of soluble subunits throughout thecytoplasm. To the best of our knowledge, how these different processes are coupled andregulated is not yet known.In a previous work we developed methods to examine and quantify the keratintransport and turnover in epithelial cells [21]; the spatial distribution of the assembledkeratin 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 toderive, from first principles, experimentally-driven mathematical models that couldyield hypothetical predictions testable in laboratories. Our approach is unique, ittranslates experimental observations and data into a series of alternative plausiblemathematical models or scenarios to further advance our understanding of the criticalparameters in keratin cycling. Fundamental to this approach is the fact that optimalparameter values for each scenario are established and out of this set, a single scenariois identified that best fits experimental observations and data.Hence, a collection of mathematical models resulting from different assumptions ofthe keratin transport, assembly and disassembly is designed to investigate the effects ofthe interplay between turnover and transport on the keratin organization. Thecollection is built as a well-designed scientific experiment by considering control andknockout of processes. To highlight and confirm the importance or existence of a givenprocess, scenarios in which the process is absent are also considered and tested. Model2/29esponses are then compared to experimental observations and data published in [21] toidentify optimal parameter values that yield the best fit of each of the models to theexperimental data. Finally, we identify, using an information-theoretic approach, thebest scenario or model given the data and candidate models under study.By employing this reductionist phenomenological approach and systematicevaluation of the different scenarios we not only confirm the proposed transport featuresof the keratin cycle and the restricted disassembly in the perinuclear region but alsofind that the assembly throughout the cytoplasm fit best to the experimental data.Furthermore, our particular approach allows us to demonstrate that the inward motionexperimentally observed is not an emergent behavior but is an inherent property of thekeratin material organization and it is due to an active transport thereby confirmingrecent experimental observations [17, 33]. Methods
Experimental data
In Moch et al. [21], the spatial distribution of the assembled keratin material inepithelial cells is measured for 15 minutes at 24 hours and 48 hours after seeding. Theshape of each epithelial cell is normalized to fit a circle of fixed radius [22]. The averagenormalized spatial distribution is calculated over 50 cells at 24 hours (resp. 84 cells at48 hours). The average speed and direction of the motion of the assembled keratinmaterial 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 toregions of assembly as
Sources and regions of disassembly as
Sinks . More details on theexperimental data can be found in [10, 21].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 thecenter of the cell and is of length 2 L with L = 22 . µm . Moch et al. recorded thefluorescence intensity of fluorescent protein-labeled keratins in cells. Assuming aproportionality between fluorescences and concentrations, and knowing from [8] themean concentration for keratin in keratinocytes, we convert fluorescence intensities toconcentrations ( µM ) as follows: Concentration = Fluorescence × (Mean Concentration/ Mean Fluorescence). In Fig. 2, the average spatial distribution, the speed of theassembled keratin, regions of assembly (Sources) and regions of disassembly (Sinks) onthe one-dimensional cross-section domain are displayed. Mathematical models
To study its organization in cells, the keratin material is categorized into a solublepool composed of the soluble keratin, and an insoluble pool representing theassembled keratin material. The state variables used to represent the soluble andinsoluble pools are: • S ( x, t ) denotes the concentration of the soluble keratin material at position x attime t , • I ( x, t ) denotes the concentration of the assembled keratin material at position x at time t .From here onwards, the one-dimensional spatial domain representing the cell is definedby Ω = { x : − L ≤ x ≤ L } ,
20 −10 0 10 200200400600800 Location ( µ m) C on ce n t r a t i on ( µ M ) Data = 24 hoursFit = 24 hours, P(x)Data = 48 hoursFit = 48 hours, f final (x) −20 −10 0 10 2001234 x 10 −3 Location ( µ m) S p ee d ( µ m / s ) Data −20 −10 0 10 20Location ( µ m) Net assembly/disassemblySourceSink
Figure 2. Experimental data adapted from [12]. µM in [15]. Circles represent raw experimentaldata; curves are fits of experimental data. Details on P ( x ) and f final ( x ) can be found in the SupportingInformation SI5. 2.2: The speed of the assembled keratin material at 24 hours. 2.3: Regions of assemblyand disassembly denoted Sources and Sinks respectively. In all figures, the cell is represented by aone-dimensional cross-section domain centered at the center of the cell; plasma membrane positions are at ± L = ± . µm , the nuclear envelope is located at ± . µm and the center of the cell is located at zero. Figure 2. Experimental data adapted from [21]. µM in [8]. Circles represent raw experimental data; curves are fits of experimentaldata. Details on P ( x ) and f final ( x ) can be found in Appendix 4. 2.2: The speed of theassembled keratin material at 24 hours. 2.3: Regions of assembly and disassemblydenoted Sources and Sinks respectively. In all figures, the cell is represented by aone-dimensional cross-section domain centered at the center of the cell; plasmamembrane positions are at ± L = ± . µm , the nuclear envelope is located at ± . µm and the center of the cell is located at zero.where x = 0 is the center of the cell and x = ± L are the boundary positions of theplasma membrane (Table 1). The general framework of the model, derived from firstprinciples based on experimental observations, takes into account the turnover andtransport for both soluble and insoluble pools and can be stated verbally as: (cid:40) Rate of change of S = Transport of S − Assembly of S in I + Disassembly of I in S, Rate of change of I = Transport of I + Assembly of S in I − Disassembly of I in S. Hence, the generalized model has the following expression; stated mathematically as: ∂S∂t = T S ( S ) − A ( S ) + D ( I ) ,∂I∂t = T I ( I ) + A ( S ) − D ( I ) , (1)where T S ( · ) (resp. T I ( · )) is the functional term describing the transport of the solublepool (resp. insoluble pool). For example, T S ( S ) = − ∂∂x J S ( x, t ) (cid:16) resp. T I ( I ) = − ∂∂x J I ( x, t ) (cid:17) where J S ( x, t ) (cid:16) resp. J I ( x, t ) (cid:17) describes the flux of the solublepool (resp. insoluble pool) at position x from the left ( x = − L ) to the right ( x = L ) attime t . The function A ( S ) (resp. D ( I )) is the assembly term (resp. disassembly term).To investigate the interplay between turnover and transport on the organization of thekeratin material several assumptions are proposed for each of the functionals T S , T I , A and D . 4/29 able 1. Model parametersParameter Definition Value (Unit) Reference L Half-length of spatial do-main Ω 22 . µm ) [21] f ( x ) Initial distribution of as-sembled keratin in cells at24 hours Eq.(17) ( µM ) [21] D S Diffusion coefficient of sol-uble pool 0 . ± .
08 ( µm /s ) [16] D I Diffusion coefficient of in-soluble pool ≤ − D S ( µm /s ) u Speed of insoluble pool 0 . − .
008 ( µm/s ) [15, 21, 34, 35] v ( x ) Almost constant speed ofactive transport of insolu-ble pool Eq.(12) ( µm/s ) [21, 34] u ( x ) Variable speed of activetransport of insoluble pool Eq.(13) ( µm/s ) [21] k ass Rate of assembly of solublepool TBD (linear model: s − ,nonlinear model: µM.s − ) k ass ( x ) Localized rate of assemblyof soluble pool Eq.(14) with k ass TBD(linear model: s − , nonlin-ear model: µM.s − ) k S Concentration for half-saturation of assembly rate(nonlinear model) TBD ( µM ) k dis Rate of disassembly of in-soluble pool into solublepool TBD (linear model: s − ,nonlinear model: µM.s − ) k dis ( x ) Localized rate of disassem-bly of insoluble into solublepool Eq.(15) or (16) with k dis TBD (linear model: s − ,nonlinear model: µM.s − ) k I Concentration for half-saturation for disassemblyrate (nonlinear model) TBD ( µM )The model parameters used in all the scenarios. (TBD = To Be Determined by fittingmodel solutions to experimental data). Modes of transport
Molecules of the soluble pool are assumed to be subjected tothe Brownian motion; the soluble pool is diffusible with a diffusion coefficient D S . The5/29unctional term of transport for the soluble pool is given by T S ( S ) := D S ∂ S∂x . (2)Passive transport for both pools is assumed in all models. Diffusion is assumed for theinsoluble pool to describe the wiggling motion of the keratin filaments in cells. Thediffusion coefficient of the insoluble pool is set to be much smaller than that of thesoluble pool: 0 < D I ≤ − D S . It is assumed that only the insoluble pool can bedriven by an active transport. Experimental evidence show that the assembledintermediate filament proteins in the form of filament precursors, squiggles or filamentsmove along microtubules and actin filaments by interacting via molecularmotors [9, 27, 33, 34, 36]. An active transport (also called the inward drift) of theinsoluble pool from the plasma membrane to the center of the cell is hypothesized basedon reports of the motion of the assembled keratin material mostly towards the nucleusin epithelial cells [15, 33, 34]. The speed of the active transport is set to be almostconstant v ( x ) everywhere or variable u ( x ) everywhere (Fig. 3). Based on experimentalobservations, both speeds are assumed to decay around the nucleus towards the centerof the cell. The variable speed u ( x ) is estimated using the profile of the average speedsmeasured in [21] (Fig. 2.2). The magnitude of the almost constant speed v ( x ) is set tobe the average value of the variable speed over the cell. Details on the derivation of theestimates of v ( x ) and u ( x ) are given in Appendix 1. Hence, the functional term for thetransport of the insoluble pool can take three different forms: T I ( I ) := D I ∂ I∂x , no drift, D I ∂ I∂x + sgn ( x ) v ( x ) ∂I∂x , inward drift with almost constant speed, D I ∂ I∂x + sgn ( x ) u ( x ) ∂I∂x , inward drift with variable speed, (3)where v ( x ) is given in (12) and u ( x ) in (13); the two functions representing the speedsof the active transport are graphed in Fig. 3. The function sgn ( x ) defined by sgn ( x ) = x > , − x < , describes the inward direction of the active transport at any location of the spatialdomain Ω centered at zero.Combining the modes of transport for the soluble and insoluble pools, three types oftransport are considered for the keratin material in the epithelial cell. Expressions of assembly / disassembly reactions
In the present work, theturnover is composed of two reactions: the assembly / aggregation / polymerization ofunits 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 thatthe assembly process is a function of the soluble pool only, whereas the disassemblyprocess is a function of the insoluble pool only. The simplest case is to consider a linearmodel that assumes linear exchanges between the two pools. The linear model can bestated as follows: ± (cid:16) A ( S ) − D ( I ) (cid:17) := ± (cid:16) k ass ( · ) S − k dis ( · ) I (cid:17) , (4)6/29
20 −10 0 10 2000.511.522.533.5 x 10 −3 Location ( µ m) S p ee d ( µ m / s ) v(x)u(x) Figure 3. Profiles of the active transport speeds u ( x ) and v ( x ) for theinsoluble pool both of which are derived from the experimental measurements in [21].The function u ( x ) is derived from the profile of speeds and v ( x ) is approximated as theaverage value of these speeds. Details on the derivation of u ( x ) and v ( x ) are given inAppendix 1.where k ass ( · ) (resp. k dis ( · )) is the rate of assembly of the soluble pool (resp.disassembly of the insoluble pool). Both rates can either be constant, k ass and k dis , orspace-dependent k ass ( x ) and k dis ( x ).In the second case, the turnover is assumed to depend on the enzymaticactivities [29]. For instance, the solubilization of the insoluble pool into soluble proteins(disassembly) is triggered by a kinase activity and the assembly of insoluble pool isinduced by the dephosphorylation of soluble proteins by a phosphatase [12]. Theturnover term is assumed to be of the Michaelis-Menten form stated as: ± (cid:16) A ( S ) − D ( I ) (cid:17) := ± (cid:18) k ass ( · ) Sk S + S − k dis ( · ) Ik I + I (cid:19) , (5)where k ass ( · ) (resp. k dis ( · )) is the maximum rate of assembly of the soluble pool (resp.disassembly of insoluble pool) and k S (resp. k I ) is the concentration at which theassembly (resp. disassembly) rate is half of k ass ( · ) (resp. k dis ( · )). WhenMichealis-Menten dynamics are used, the model is called nonlinear. Both rates caneither be constant or space-dependent describing the intracellular localization of thepost-translational modification enzymes. Profiles of assembly and disassembly rates
As previously mentioned, theassembly and disassembly rates, k ass ( · ) and k dis ( · ), used in the linear and nonlinearmodels can be constant or space-dependent functions. The profile (shape) of thespace-dependent function k ass ( x ) is derived from the spatial profile of regions ofassembly (Sources) measured in [21] (Fig. 2.3). Details of the derivation of the Sources k ass ( x ) are given in Appendix 2.Two types of shapes for k dis ( x ) are assumed to represent two types of localization ofthe disassembly in cells. First, similarly to the assembly rate, the profile of thedisassembly rate is deduced from the experimental data published in [21] (Fig. 2.3); thespatial profile of the disassembly regions, Sinks, is used to build the shape of the firstspace-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-functionis designed to describe this assumption. This second space-dependent disassembly rateis called of type Mollify. Details on the derivation of the two k dis ( x ) rates are given inAppendix 3. 7/29or the sake of illustration, the two profiles of k ass ( · ) and the three profiles of k dis ( · )are shown in Fig. 4. −20 −10 0 10 200123456 x 10 −3 A sse m b l y r a t e Location ( µ m) k ass k ass (x): Source −20 −10 0 10 200
234 x 10 −3 D i sasse m b l y r a t e Location ( µ m) k dis k dis (x): Sinkk dis (x): Mollify Figure 4. Possible profiles for the assembly and disassembly rates. k ass = 10 − s − whenthe linear model is used (resp. µM.s − when the nonlinear model is used) and k ass ( x ) is computed using(14). 4.2: k dis = 10 − s − when the linear model is used (resp. µM.s − when the nonlinear model isused) and k dis ( x ) are computed using (15) for Sinks and (16) for Mollify, respectively. Details on thederivation of k ass ( x ) and k dis ( x ) rates are given in the Supporting Information SI3 and SI4, respectively. Figure 4. Possible profiles for the assembly and disassembly rates. k ass = 10 − s − when the linear model is used (resp. µM.s − when the nonlinear modelis used) and k ass ( x ) is computed using (14). 4.2: k dis = 10 − s − when the linear modelis used (resp. µM.s − when the nonlinear model is used) and k dis ( x ) are computedusing (15) for Sinks and (16) for Mollify, respectively. Details on the derivation of k ass ( x ) and k dis ( 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 givethe same total amount of assembly / disassembly over the spatial domain Ω.Accounting for the three modes of transport and considering the turnover describedby either linear and nonlinear models with 6 possible combinations of the profiles for theassembly and disassembly rates, a collection of 36 scenarios (models) is defined; eachscenario follows the general form stated by system (1). Details of the 36 scenarios aregiven in Fig. 5. All scenarios are considered with the same initial conditions given by S ( x, t ) = 0 . . f ( x ) ,I ( x, t ) = f ( x ) , for all x ∈ Ω , (6)where f ( x ), defined in (17), is chosen to be a mollified version of the profile of theassembled keratin at time t = 24 hours averaged over all the normalized cells (Fig. 6).Details of the derivation of f ( x ) are given in Appendix 4. Initial conditions describethe observed fact that the soluble pool (resp. insoluble pool) represents 5% (resp. 95%)of the total keratin material [4]. All scenarios are considered with boundary conditionsdescribing the impermeability of the plasma membrane for the keratin material J S ( ± L, t ) = 0 = J I ( ± L, t ) , t ≥ t , (7)where J S is the flux of the soluble pool and J I is the flux of the insoluble pool as definedbelow system (1). The model parameters used in all the scenarios are listed in Table 1.8/29 cenario i Diffusion No drift Cst drift Var. drift Lin. Nonlin. Cst Ass. Sources Cst Diss. Sinks Mollify
Figure 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 . Figure 5. The 36 scenarios to be considered.
Top: Number in parentheses is thescenario index i . Bottom: The numerical value 1 denotes that the process of interest isin Scenario i . Comparison between mathematical models and experimentaldata
Parameter estimation
Let p denote the set of all the model parameters for eachscenario. To estimate the optimal set of parameter values p for each scenario, the 9/29
20 −10 0 10 200100200300400500600700 Location ( µ m) C on ce n t r a t i on ( µ M ) f (x) Figure 6. Initial profile of the insoluble pool f ( x ) defined in (17).solution of each scenario is compared to experimental data using the following objectivefunction (an estimation of the distance between the experimental data and the modelresponse): Ψ( p ) = N x (cid:88) i =1 (cid:104) I ( x i , t final , p ) − f final ( x i ) (cid:105) , (8)where t final is equal to 48 hours. The experimental data at t final is represented by f final ( x ) and the model solution at t final of the considered scenario evaluated with theparameter values p is I ( x, t final , p ). To obtain the model solution for each scenario, thecorresponding system is numerically integrated from t = 24 hours to 50 hours using theMATLAB solver for partial differential equations, pdepe [20]. Solutions are computed at N x locations of the spatial domain Ω with N x = 200. To carry out computations, theraw data in which the concentration of the assembled keratin material is known at 623spatial locations is approximated by f final ( x ) defined in Appendix 4 by (18); f final ( x )is the fit of the average profile of the assembled keratin measured after 48 hours on thenormalized cells (see the black profile in Fig. 2.1).The estimation of the parameter values for the 36 scenarios, i.e. the determination ofparameter values ˆ p that provide the best fit to the experimental data, is done byminimizing the objective function, Ψ( p ), such thatΨ(ˆ p ) = min p Ψ( p ) . (9)Minimization of the objective function is done by a parallelizable genetic algorithmdescribed in [7].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 ass and k dis are estimated. When the nonlinear model is considered, four parameters, k ass , k S , k dis and k I , are estimated for each scenario. In all scenarios, the diffusion coefficients for thesoluble and insoluble pools as well as the active transport speeds v ( x ) and u ( x ) areconsidered as fixed/known parameters or functions and are set to values measured inprevious studies [16, 21, 33] (Table 1). The diffusion coefficients are taken as D S = 0 . µm /s and D I = 9 . × − D S , and the “constant speed” u in v ( x ) asdefined in (12) is set to be u = 0 . µm/s [16, 34]. Model selection
In order to select the best model out of the 36 scenarios, we use aninformation-theoretic approach. Specifically, we employ the Akaike information criterion10/29AIC) [3, 13] to select from all the scenarios, the best model that best capturesexperimental observations given the collection of models considered in this study. Sincewe use the Least Squares principle to estimate the values of the constant freeparameters (2 free parameters for linear models and 4 for nonlinear models), the Akaikecriterion for Scenario i , AIC i , takes the following form: AIC i = N x ln (cid:0) ˆ σ (cid:1) + 2 K, (10)where ˆ σ = Ψ i (ˆ p ) N x is the estimate of the variance with Ψ i (ˆ p ) the residual (9) estimatedfor Scenario i and N x the size of the sample ( N x =200). K is the number of estimatedparameters; i.e. the number of free parameters in Scenario i plus one for the estimate ofthe variance ( K = 3 for linear model and K = 5 for nonlinear models). The scenariowith the lowest AIC value is the best model. The AIC selects a model with the leastnumber of parameters that best fits experimental observations. To rank and comparescenarios the Akaike weights w i are calculated and these are known as the weights ofevidence in favor of Scenario i being the actual best model given the experimental dataand the collection of scenarios considered. The Akaike weights are expressed as follows: w i = exp ( − / i ) (cid:80) Rr =1 exp ( − / r ) , with ∆ i = AIC i − min i AIC i (11)where R represents the number of scenarios considered ( R = 36) and ∆ i being thedifference in AIC with respect to the AIC of the preferred scenario min i AIC i = AIC min .It must be noted that the Akaike weights w i sum to 1 and are interpreted as theprobability that Scenario i is the best model given the experimental data and thecollection of scenarios considered. The models, ranked from the largest to the smallestAkaike weights, whose Akaike weights sum to 0.95 form the confidence set of the modelsthat captures, more faithfully, experimental data. The ratio of the Akaike weights w i /w j (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 ofall scenarios involving the process of interest. We will denote by w + ∗ the sum of theAkaike weights of the scenarios including the process of type ∗ . This sum can beinterpreted as the probability that the process of type ∗ is the best type of the processgiven 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 foreach scenario is found by minimizing the objective function (9). Secondly, the Akaikeinformation criterion (10) is used to select the best of the best fits out of all thescenarios.The best fit for each of the 36 scenarios is presented in Fig. 7 in the following order: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 turnoverterms. • Rows 1 and 4 are scenarios with constant disassembly rate. • Rows 2 and 5 are scenarios with localized disassembly rate of type Sinks.11/29
Rows 3 and 6 have scenarios with localized disassembly rate of type Mollify.2. Column-wise: • The first two columns (1 and 2) include scenarios with diffusion alonewithout 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 typeSources.According to the
AIC i values, the best model, the best of the best fits, is identifiedas being Scenario 21, AIC = min i AIC i (3 rd column of Table 2). Since Scenario 21’sAkaike weight is over 0 .
95 (5 th column of Table 2), it is the only model out of the setconsidered that satisfies the confidence criterion; Scenario 21 matches very wellexperimental observations and data. The second best model is Scenario 31. Theevidence ratio of scenarios 21 and 31 w /w is equal to 36 . × ; i.e. Scenario 21 isabout 36 millions times more likely than Scenario 31 to be the best model given theexperimental data and the collection of models considered. We consider this to bestrong evidence in support of Scenario 21. 12/29 able 2. Results of the model selection for the best fit of each of the 36scenarios.Scenario i K AIC i ∆ i w i Rank . × −
232 5 2339.200 836.5561 0 3033 5 2119.447 616.8033 0 2434 5 2074.025 571.3811 0 2235 5 1834.186 331.5420 0 636 5 2143.595 640.9505 0 25The numerical value 0 in the 5 th column denotes a value lower than 10 − . In the 6 th column, “Rank” denotes the ranking of scenarios, i.e. the descending order of Akaikeweights w i . The use of weights allows the ”quantitative” comparison of the adequacy ofscenarios because of the normalization to 1.Scenario 21 is characterized by an inward drift with an almost constant speed for theinsoluble pool, turnover terms of Michaelis-Menten type, a constant assembly rate and adisassembly rate of type Mollify. The profiles of assembly and disassembly rates aredisplayed in Fig. 8. For this scenario, the estimated optimal parameter values are k ass = 9 . µM/s , k S = 570 . µM , k dis = 0 . µM/s leading to k max = 1 . µM/s in (16) and k I = 976 . µM . Based on the Michaelis-Menten constants for the assembly13/29nd disassembly processes, since k S < k I , an enzyme that would be involved in thesolubilization of the assembled keratin material requires a higher substrateconcentration to achieve a given reaction speed than an enzyme that would be involvedin the assembly of soluble proteins. Snapshots of the soluble and insoluble pool profilestaken every 30 minutes from t to t final are displayed in Fig. 8. After only 2 hours, thesolution stabilizes to its final profile. Scenario 21 preserves the repartition of the keratinmaterial between the soluble and insoluble pools over time, about 95% of the keratinmaterial is assembled to form the insoluble pool. Finally, the characteristic time scalesof the passive transport ( τ Diffusion ), active transport ( τ Drift ) and turnover ( τ Reaction )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 timescales of the processes included in Scenario 21 are ordered as follows: τ Reaction (10 s ) < τ SDiffusion (10 s ) < τ Drift (10 s ) < τ IDiffusion (10 s ) . Using Peclet’s number (
P e = τ Diffusion /τ Drift (cid:29) Da = τ Drift /τ Reaction (cid:29) ω i of all the scenarios is given inthe 6 th column of Table 2. It is worthwhile remarking that differences in AIC i , ∆ i ,could have been used for ranking purposes; in this case, we would have obtained thesame ranking. The rationale of using ω i is that it allows us to quantify how preferablyeach candidate model is via the normalization to 1. Importance of the type of process
Using the sums of the Akaike weights, w + ∗ , we investigate several questions to evaluatethe 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 ofprocesses considered. For instance, considering the transport process, the collection ofthe 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 from13 to 24) and scenarios with a variable speed drift (Scenarios from 25 to 36). The sumof the Akaike weights of each category is then calculated, compared and ordered tocharacterize which type is more likely to be present. In what follows the sign > denotesthe 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 w + NoDrift w + CstDrift w + V arDrift w + Linear w + Nonlinear w + CstAss w + Source w + CstDis w + Sink w + Moll =0 =1 = 2 . × − = 0 =1 =1 =0 = 2 . × − =0 =1 w + ∗ 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 withvariable 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 typeSinks vs Localized of type Mollify. 14/29 able 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. w + CstAssDis w + CstAssSink w + CstAssMoll w + SourceDis w + SourceSink w + SourceMoll = 2 . × − =0 =1 =0 =0 =0 Combinations of assembly/disassembly rate profiles - Constant assembly and disassembly rates vs constantassembly rate and disassembly rate of type Sinks vs constant assembly rate and disassembly rate of typeMollify vs assembly rate of type Sources and constant disassembly rate vs assembly rate of type Sources anddisassembly 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: w + ∗ 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 incells 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), driftwith an almost constant speed for the insoluble pool or drift with a variable speedfor the insoluble pool ? Each category includes 12 scenarios (see 3 rd to 5 th columnsin Fig. 5). From the Akaike weights, it is obvious that scenarios includingdiffusion only (with no drift) are not supported by experimental data since w + NoDrift = 0 (Table 3). Diffusion alone is not enough to explain the experimentaldata. Given the data and the model collection, the type of transport can beordered as follows (Table 3):
Drift with almost constant speed > Drift with variable speed (cid:29)
No drift .15/29
What is the most likely type of turnover between the soluble andinsoluble pools given the experimental data and the model collectionconsidered?
Is linear exchange or Michaelis-Menten type (underlying an enzymeactivity) more likely ? Each category includes 18 scenarios (see 6 th to 7 th columnsin Fig. 5). From Table 3, enzymatic activities are more likely to be present. Hence Nonlinear exchange (cid:29)
Linear exchange. • What is the most likely rate profile of the assembly process of thekeratin material in cells given the experimental data and the modelcollection considered?
Is a constant assembly rate or assembly mainly localizedat the cell membrane periphery more likely ? Each category includes 18 scenarios(see 8 th to 9 th columns in Fig. 5). From Table 3, the non-localized rate of assemblyis preferred to the rate of the localized assembly; w + CstAss > w + Source . Hence
Non-localized assembly (cid:29)
Assembly localized at cell membrane periphery. • What is the most likely rate profile of the disassembly process of thekeratin material in cells given the experimental data and the modelcollection considered?
Is a constant disassembly, disassembly of type Sinks ordisassembly localized around the nucleus more likely ? Each category is composedof 12 scenarios (see 10 th to 12 th columns in Fig. 5). From Table 3, none of thescenarios that include the disassembly rate of type Sinks is supported byexperimental data. Disassembly localized around the nucleus > Non-localized disassembly (cid:29)
Sinks-type profile. • What is the most likely combination of the assembly and disassemblyrate profiles given the experimental data and the model collectionconsidered?
Is a combination of constant assembly and disassembly rates, aconstant assembly rate and disassembly rate of type Sinks, a constant assemblyrate and disassembly rate of type Mollify, an assembly rate of type Sources andconstant disassembly rate, an assembly rate of type Sources and disassembly rateof type Sinks or an assembly rate of type Sources and disassembly rate of typeMollify more likely?
Now, interactions between the assembly and disassemblyprocesses are investigated. Six categories are defined. In each category, there are 6scenarios as listed in the top of Table 4. Overall, as in Scenario 21, thecombination of non-localized assembly and disassembly localized around thenucleus is preferred (Table 4).
Non-localized assembly and disassembly localized around the nucleus > Non-localized assembly and disassembly (cid:29)
Any other combination.
Discussion
The primary aim of the present work is to investigate, through mathematical modelingdriven by experimental observations, mechanisms contributing to the organization ofthe keratin material in epithelial cells and subsequently to test the biological model forthe keratin dynamics proposed by Leube and Windoffer et al. in 2011 [17, 33]. Fromfirst principles, we formulate a collection of mathematical models capturing variouscombinations of biological processes describing the spatio-temporal dynamics of the16/29eratin material in epithelial cells. By using techniques in parameter estimation, we findoptimal reaction kinetic parameter values for each model that best capturesexperimental observations. We go one step further and employ an information-theoreticapproach for model selection to determine the best model among all the best fit modelsthat captures the key processes of the experimental data.As previously highlighted, the model framework and the description of processesconsidered are driven by experimental data and observations. Experimental data used inthis study provide the spatial distribution of the assembled keratin concentration at 24and 48 hours. As only macroscopic information on the keratin organization is available,the models only describe the keratin material as soluble or assembled and considerexchanges between these soluble and insoluble pools combined with transport events.For instance, as experimental data identify the existence of regions with preferentialassembly and disassembly (Fig. 2.3), the assembly and/or disassembly processes areassumed to be localized in some scenarios. When observing the perpetual inward motionof the keratin network (see S1 Video), the “natural assumption” for the transport of theassembled keratin would be the existence of an inward active transport. However, wedecide 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 patternformation (see some examples in [19]). Secondly, to reinforce the predictive power of ourstudy. In our approach, we do not only design mathematical models but also calibratemodels (parameter estimation) and then evaluate (model selection) how each modelperforms. 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. Ourstrategy of considering “unrealistic scenarios” and “realistic scenarios” in the collectionof models to be evaluated mimics the biological experimental protocols such as knockoutand controls. For instance, it has been shown that the keratin assembly / disassemblydepend on post-translational modifications of keratins due to enzymatic activities [12].In our study, a total of 36 scenarios that divide into 18 scenarios with a turnover of typelinear (no enzymatic activities) and 18 scenarios with a turnover of type nonlinear(enzymatic activities) are investigated. There is a one-to-one correspondence betweenthe 18 scenarios of the 2 groups. When model performances (Akaike weights) arecompared, scenarios having a nonlinear turnover perform better than the correspondingones with a linear turnover. Hence, the enzymatic activity is detected by the modelselection as existing but also preferable. This outcome allows us to judge that ourapproach 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 bestfits. Scenario 21 is in good agreement with the biological model proposed in [17, 33].Scenario 21 and the biological model share the following common key features: diffusionof the soluble pool, inward active transport of the insoluble pool and disassembly of theinsoluble pool localized around the nucleus . Moreover, as experimentally observed,Scenario 21 preserves the repartition of the keratin material over time; the keratinmaterial is mainly insoluble in epithelial cells. Both the proposed biological model andScenario 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 12in Table 2, 1 st and 2 nd columns in Fig. 7 and 1 st to 3 rd columns in Table 3), we go onestep further and show that the active transport is a requirement to explain theexperimental data and that the experimentally observed inward motion of theassembled 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 thekeratin dynamics are mainly controlled by the active transport in Scenario 21. 17/29nterestingly, Scenario 21 concurs with the proposed biological model about theexistence of a disassembly localized in the perinuclear region. An experimental protocolmust now be developed to further work out the details of this conclusion. It is worthpointing out that some characteristics of the proposed biological model were not testeddue to the form of the mathematical models. For instance, the scenarios considered heredescribe the turnover in terms of the assembly of soluble proteins and disassembly ofaggregated proteins. In the biological model, the nucleation of filaments, i.e. theinitiation 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 wasincluded 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 nucleationterm was subsequently dropped. Investigation, by mathematical modeling, of thenucleation 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 th and6 th columns in Fig. 7), the preferred combination of assembly and disassembly rateprofiles is a constant rate of assembly and disassembly as in Scenario 31, which is thesecond best model. The best model, Scenario 21, and the second best model, Scenario31, share the non-compartmentalization of the assembly process. However, in Scenario21 the active transport has an almost constant speed whereas in Scenario 31 the speedis variable. For the disassembly, this process is non-localized in Scenario 31 whereas it islocalized 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 ∈ [ − . , .
5] in Fig. 3) “compensates” for the non-compartmentalization of thedisassembly process. After comparing the features of Scenarios 21 and 31, we inspecttheir profiles. The profile obtained with Scenario 21 fits very well the experimental dataon non-perinuclear locations (see for x / ∈ [ − . , .
5] in Fig. 7.33). On the other hand,the profile resulting from Scenario 31 fits very well experimental data on perinuclearlocations (see for x ∈ [ − . , .
5] in Fig. 7.23). As an improvement of Scenario 21, wewould 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, alocalized assembly rate having the spatial profile deduced from Sources and a localizeddisassembly rate whose shape is derived from the profile of the Sinks. All thecharacteristics extracted from the experimental data are included in Scenario 29 (linearmodel) and Scenario 35 (nonlinear model). If one wanted the best model that takes intoaccount all the experimental features, then Scenario 29 or 35 would be the best scenario.Scenario 29 is ranked 27 th and Scenario 35 is ranked 6 th (Table 2). Similarly to thegeneral trend, for the same set of assumptions, the use of the Michaelis-Menten typeturnover term generally provides a better representation of experimental data than theuse of the linear turnover term. The evidence ratio of Scenarios 35 and 29 ω /ω isridiculously large; Scenario 35 is much more adequate to represent the experimentaldata than Scenario 29. However, Scenario 35 is still only ranked 6 th . Thefailure/mismatch of Scenario 35 might be explained by the redundancy of theinformation existing in the variable speed and the net assembly/disassembly regionprofiles extracted from the experimental data. Furthermore, according to the Akaikeweights, disassembly rates of type Sinks (the type deduced from experimental data) areless 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 aconsequence of our too conservative interpretation of the regions of preferentialassembly or disassembly. In our assumptions, in zones of preferential assembly, thedisassembly rate is set to be about null and vice versa (Fig. 2.3). 18/29n summary, to model the keratin dynamics in epithelial cells we characterized thekeratin material into two pools, the soluble and the insoluble pool; and, the eventsconsidered are the turnover and transport for both pools. The modeling assumptionsused, for instance, the diffusion of the soluble pool, are based on biological observationsand experimental data. The collection of the models considered in this study is designedto answer a set of questions such as “what is the mode of transport of the keratinmaterial in epithelial cells?”. After optimizing parameter values, model selection andevaluation methods applicable to non-nested models are used to discriminate betweenthe candidate models, identify the best model and quantify how models underconsideration are adequate to explain the experimental data. Note that the ranking ofthe 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 thecandidate models considered here. For instance, considering other hypotheses such asthe non-negligence of the anterograde motion of the assembled keratin material or thestabilization (protection against disassembly) of the keratin filaments involved in thenuclear cage would have led to a different collection of scenarios in which our bestscenario could have failed to be the best one. Furthermore, as some modelingassumptions are directly derived from experimental data, missing information inexperimental data might have been prejudicial for the correct approximation / modelingof processes; for instance, missing information about the speed of assembled keratin inthe cell periphery result in small values for the variable speed close to the plasmamembrane (Fig. 2.2). Keeping in mind the limitations of our approach, importantconclusions have been reached such as • an active transport of the assembled keratin material is required, therebyconfirming recent experimental observations [17, 33], • enzymatic activities regulating the assembly / disassembly are more likely tooccur, • 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 isidentified, this scenario is in good agreement with the biological model previouslyproposed in [17, 33]. Interestingly, the best scenario supports the perinuclearlocalization of the disassembly process hypothesized in the biological model.Mathematical models of the keratin intermediate filament organization in cells werepreviously proposed [1, 14, 24, 25, 31]; however, none of these models included the effectsof transport of the assembled material on its organization. More importantly, only thebehavior of models were validated qualitatively, no comparisons to experimental datawere carried out. It is worth noting that other studies of the intermediate filamentsdynamics combining mathematical modeling and experimental data approaches werecarried out but on neurofilaments in neurons (see for example [2, 18]). Conclusion
Given the experimental data published in [21], through modeling and simulations, weinvestigate the effects of the interplay between turnover and transport on the keratinspatio-temporal organization in epithelial cells. Out of all the scenarios investigated, ascenario strongly supported by experimental data is found that best captures most ofthe hallmarks of the experimental observations. This scenario predicts the diffusion ofsoluble keratin, an inward active transport of the assembled keratin and the disassemblylocalized around the nucleus triggered by enzymatic activities as well as the assembly19/29rocess that is non-compartmentalized over the cell. The value of our models isreflected in their predictive nature; first, the approach predicts the localized disassemblyat the perinuclear region and second, that the experimentally observed inward motion isnot an emergent behavior but that it is an inherent property of the organization of thekeratin 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 isassumed to decay to almost zero around the nucleus. The nuclear envelope positions areat x = ± . µm . In order to describe the decay of the speed around the nucleus suchthat it is almost constant, the following function is used: v ( x ) = u (1 − exp ( − ax )) , (12)with a = 0 .
05 and u being in the range given in Table 1. For numerical simulations, u isset to 0 . µm/s which represents the average value of speeds measured in [21].For the variable speed case, a symmetrical function over the spatial domain Ω issought to describe the speed u ( x ). Averaging over the symmetrical spatial locationsvalues of the average speed measured in [21] (Fig. 2.2) and curve fitting with a sum ofGaussian functions of these values, an estimate of u ( x ) is obtained and expressed asfollows: u ( x ) = (cid:88) i =1 a i × exp ( − (( x − b i ) /c i ) ) , (13)with the following coefficients: a = 0 . b = 17 . c = 7 . a = 0 . b = − .
41 and c = 7 . u ( x ) is almost zero at the nucleus-locationsand is symmetrical around the center of the cell (Fig. 9).For numerical simulations, when needed, the following function s ( x ) = 21 + e − x/a − a = 1 is used as a “smooth analogue” of sgn ( x ). Appendix 2
To obtain the space-dependent function k ass ( x ) the profile of the regions of assembly(Sources) published in [21] (Fig. 2.3) is first made symmetrical by averaging values atthe symmetrical spatial locations and then the symmetrical Sources profile is fittedusing the sum of Gaussian functions: k ass ( x ) = k max (cid:88) i =1 a i × exp ( − (( x − b i ) /c i ) ) + k baseline . (14)The best fit is obtained with the following coefficients: a = 1 . × , b = 20 . c = 2 . a = 8 . × , b = − . c = 1 . a = 1 . × , b = − . c = 1 . k max determines the maximal value of thepeaks. When k max = 2 L ( k ass − k baseline ) √ π (cid:80) i =1 a i | c i | with k baseline = k ass × − , the totalamount of assembly over the cell is the same as that of the case of a constant assemblyof level k ass . The value k ass in k max is determined by fitting model solutions toexperimental data. An illustration of the shape of k ass ( x ) is given in Fig. 10. 20/29 ppendix 3 To obtain the first space-dependent function k dis ( x ), the profile of regions ofdisassembly (Sinks) published in [21] (Fig. 2.3) is made symmetrical and fitted using thesum of Gaussian functions: k dis ( x ) = k max (cid:88) i =1 a i × exp ( − (( x − b i ) /c i ) ) . (15)The best fit of the symmetrical profile of Sinks is obtained with the followingcoefficients: a = − . × , b = − . c = 4 . a = − . × , b = − . c = 2 . a = 2 . × , b = − . c = 4 . a = − . × , b = 5 . c = 5 . a = − . × , b = − . c = 5 . a = 3 . × , b = 11 . c = 4 . a = − . × , b = 14 . c = 2 . a = − . × , b = 11 .
46, and c = 4 . k max determines the amplitude of the peaks.When k max = 2 Lk dis √ π (cid:80) i =1 a i | c i | , the total amount of disassembly over the cell is the sameas that of the case of the constant function of level k dis . The value k dis in k max isdetermined by fitting model solutions to experimental data. An illustration of the shapeof k dis ( x ) obtained from Sinks is given in Fig. 11.A second function is hypothesized to represent the localized disassembly k dis ( x )around the nucleus. A mollified piecewise function is used and is of the form k dis ( x ) = k baseline , x < a − (cid:15), − k max (cid:15) ( a − a ) x − k max ( (cid:15) − a )2 (cid:15) ( a − a ) x + − ( (cid:15) − a ) k max + 4 (cid:15) ( a − a ) k baseline (cid:15) ( a − a ) , a − (cid:15) ≤ x < a + (cid:15),k baseline + k max x − a a − a , a + (cid:15) ≤ x < a − (cid:15), − k max (cid:15) ( a − a ) x + k max ( (cid:15) + a )2 (cid:15) ( a − a ) x − (( (cid:15) − a ) + 4 (cid:15)a ) k max + 4 (cid:15) ( a − a ) k baseline (cid:15) ( a − a ) , a − (cid:15) ≤ x < a + (cid:15),k baseline + k max , a + (cid:15) ≤ x < a − (cid:15),k max (cid:15) ( a − a ) x + k max ( (cid:15) − a )2 (cid:15) ( a − a ) x + (( (cid:15) + a ) − (cid:15)a ) k max + 4 (cid:15) ( a − a ) k baseline (cid:15) ( a − a ) , a − (cid:15) ≤ x < a + (cid:15),k baseline + k max − k max x − a a − a , a + (cid:15) ≤ x < a − (cid:15),k max (cid:15) ( a − a ) x − k max ( (cid:15) + a )2 (cid:15) ( a − a ) x + ( (cid:15) + a ) k max + 4 (cid:15) ( a − a ) k baseline (cid:15) ( a − a ) , a − (cid:15) ≤ x < a + (cid:15),k baseline , a + (cid:15) ≥ x, (16)where a = − µm , a = − . µm , a = 7 . µm , a = 15 µm and (cid:15) = 1. When k max = 2( k dis − k baseline ) with k baseline = k dis × − , the total amount of disassemblyover the cell is the same as that of the case of a constant disassembly of level k dis ; k dis is determined by fitting the model solutions to experimental data. An illustration of theshape of (16) is given in Fig. 4. 21/29 ppendix 4 To define the initial condition f ( x ), the experimental profile [21] of the assembledkeratin material measured at 24 hours is used (gray circles in Fig. 2.1). The polynomial P ( x ) = (cid:88) i =0 p i x i for which p = 5 . × − , p = 3 . × − , p = − . × − , p = − . × − , p = 0 . p = 2 . × − p = 0 . p = 1 . × − and p = 506 . P ( x )is then modified to obtain a function satisfying the boundary conditions defined in (7).The initial condition f ( x ) describing the profile of the insoluble pool at 24 hours isexpressed as follows: f ( x ) = c, x < a − (cid:15),P ( a + (cid:15) ) − c (cid:15) x + ( (cid:15) − a )( P ( a + (cid:15) ) − c )2 (cid:15) x + ( (cid:15) − a ) P ( a + (cid:15) ) + (3 (cid:15) + 2 a (cid:15) − a ) c (cid:15) , a − (cid:15) ≤ x < a + (cid:15),P ( x ) , a + (cid:15) ≤ x < a − (cid:15), − P ( a − (cid:15) ) − c + 2 (cid:15)D ( a − (cid:15) )4 (cid:15) x − ( (cid:15) − a )( P ( a − (cid:15) ) − c ) − (cid:15)a D ( a − (cid:15) )2 (cid:15) x + (3 (cid:15) + 2 a (cid:15) − a ) P ( a − (cid:15) )4 (cid:15) + 2 (cid:15) ( (cid:15) − a ) D ( a − (cid:15) ) + ( (cid:15) − a ) c (cid:15) , a − (cid:15) ≤ x < a + (cid:15),c, a + (cid:15) ≤ x, (17)with a = − µm , a = − a , (cid:15) = 1, c = 50 µM and D ( x ) = (cid:88) i =1 ip i x i − . The function f ( x ) is graphed in Fig. 6.The average profile of the assembled keratin material computed on the normalizedcells after 48 hours of seeding [21] is represented by: f final ( x ) = (cid:88) i =0 p i x i , (18)where p = − . p = 2 . × − , p = 0 . p = 1 . × − and p = 604 . f final ( x ) is graphed in Fig. 2.1 (black curve). Appendix 5
The system corresponding to Scenario 21 is now nondimensionalized to estimate thecharacteristic time scale of each of the processes involved. Consider new independentvariables z and τ and new dependent variables ˜ S ( z, τ ) and ˜ I ( z, τ ) defined such that: x = Az, t = Bτ, S ( x, t ) = s ˜ S ( z, τ ) , I ( x, t ) = i ˜ I ( z, τ ) . Scenario 21 takes then the following form: ∂ ˜ S∂τ = BD S A ∂ ˜ S∂z − Bs (cid:32) k ass ˜ Sk S /s + ˜ S − k dis ( Az ) ˜ Ik I /i + ˜ I (cid:33) ,∂ ˜ I∂τ = B(cid:15)D S A ∂ ˜ I∂z + BuA sgn ( z ) f ( Az ) ∂ ˜ I∂z + Bi (cid:32) k ass ˜ Sk S /s + ˜ S − k dis ( Az ) ˜ Ik I /i + ˜ I (cid:33) , f ( Az ) being the functional component of v ( x ) (12) defined as f ( Az ) = (1 − exp ( − a ( Az ) )) with a = 0 .
05 and k dis ( Az ) being the function k dis ( · )defined in (16) and evaluated at Az . Taking A = (cid:96) ( (cid:96) = 2 L ), B = (cid:96)/u , s = k S and i = k S then x = (cid:96)z and t = (cid:96)u τ and the adimensional version of the system is ∂ ˜ S∂τ = 1 (cid:15)P e ∂ ˜ S∂z − Da (cid:32) ˜ S S − k dis ( Az ) /k ass ˜ Ik I /k S + ˜ I (cid:33) ,∂ ˜ I∂τ = 1
P e ∂ ˜ I∂z + sgn ( z ) f ( Az ) ∂ ˜ I∂z + Da (cid:32) ˜ S S − k dis ( Az ) /k ass ˜ Ik I /k S + ˜ I (cid:33) . The constant
P e = τ Diffusion τ Drift = (cid:96) /D I (cid:96)/u = u(cid:96)(cid:15)D S is the P´eclet number of the insolublepool that compares active transport and diffusion processes. The characteristic length (cid:96) is chosen to be equal to the length of the spatial domain Ω, (cid:96) = 2 L . The Damk¨ohlernumber Da = τ Drift τ Reaction compares reaction and active transport processes. Thecharacteristic active transport time is τ Drift = (cid:96)/u as defined in the P´eclet number.The characteristic reaction time is τ Reaction = 1 /κ for which κ = k ass k S . Note that sgn ( z ) f ( Az ) is about ± k ass = 9 . µM/s , k S = 570 . µM , k I = 976 . µM/s and k max = 1 . µM/s in (16)that correspond to a constant disassembly rate of level k dis = 0 . µM/s . Hence, fromthe definition of k max , k dis ( Az ) can be approximated by k dis = 0 . D S = 0 . µm /s , D I = (cid:15)D S with (cid:15) = 9 . × − and u = 0 . µm/s . The following estimates are then obtained: • Adimensional parameters: k dis ( Az ) /k ass = k dis /k ass = 0 . k I /k S = 1 . • Diffusion time scale: τ Diffusion = (cid:96) /D . – For the soluble pool: τ SDiffusion = (cid:96) /D S ≈ . × s , – For the insoluble pool: τ IDiffusion = (cid:96) / ( (cid:15)D S ) ≈ . × s ; • Active transport time scale: τ Drift = (cid:96)/u = 1 . × s (time the fluid needs toflow through the characteristic length); • Reaction time scale: k S /k ass ≈ s (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´eclet andDamk¨ohler numbers are used: • P e = ( u(cid:96) ) / ( (cid:15)D S ) ≈ (cid:29)
1. Diffusion is small compared to active transport forthe insoluble pool, the transport of the assembled keratin material by drift isfaster than by diffusion, active transport is the dominant mode of transport; • Da = ( (cid:96)k ass ) / ( uk S ) ≈ (cid:29)
1. Active transport time scale is greater than thereaction time scale; the overall process in controlled by active transport. 23/29 cknowledgments
This work was initiated in June 2013 while SP was visiting the Institute of Molecularand Cellular Biology at RWTH Aachen University (Aachen, Germany). The authorsacknowledge suggestions from anonymous referees and the editor for the critical reviewthat helped improve the quality of the manuscript.
References
1. M. Beil, S. L¨uck, F. Fleischer, S. Portet, W. Arendt, and V. Schmidt. Simulatingthe formation of keratin filament networks by a piecewise-deterministic markovprocess. J. Theor. Biol., 256:518–532, 2009.2. A. Brown, L. Wang, and P. Jung. Stochastic simulation of neurofilamenttransport in axons: the “stop-and-go” hypothesis. Mol. Biol. Cell., 16:4243–4255,2005.3. K. Burnham and D. Anderson. Model Selection and Multimodel Inference. APractical Information-Theoretic Approach. Springer Verlag, second edition, 2002.4. C.-F. Chou, C. Riopel, L. Rott, and B. Omary. A significant soluble keratinfraction in “simple” epithelial cells lack of an apparent phosphorylation andglycosylation role in keratin solubility. J. Cell Sci., 105:433–444, 1993.5. B. Chung, J. Rotty, and P. Coulombe. Networking galore: intermediate filamentsand cell migration. Curr. Opin. Cell Biol., 25:600–612, 2013.6. P. Coulombe, Kerrns M., and E. Fuchs. Epidermolysis bullosa simplex: aparadigme for disorder of tissue fragility. J. Clin. Inv., 119:1784–1793, 2009.7. R. Dorsey and W. Mayer. Genetic algorithms for estimation problems withmultiple optima, nondifferentiability and other irregular features. J. Bus. Econ.Stat., 13:53–66, 1995.8. X. Feng, H. Zhang, J. Margolick, and P. Coulombe. Keratin intracellularconcentration revisited: implications for keratin function in surface epithelial. J.Inv. Derm., 113:850–853, 2013.9. B.T. Helfand, L. Chang, and R.D. Goldman. Intermediate filaments are dynamicand motile elements of cellular architecture. J. Cell Sci., 117:133–141, 2004.10. G. Herberich, R. Windoffer, R. Leube, and T. Aach. 3d segmentation of keratinintermediate filaments in confocal laser scanning microscopy. In Engineering inMedicine and Biology Society, EMBC, 2011 Annual International Conference ofthe IEEE, pages 7751–7754, Aug 2011.11. M. Homberg and T. Magin. Beyond expectations: novel insights into epidermalkeratin function and regulation. Int. Rev. Cell Mol. Biol., 311:265–306, 2014.12. I. Izawa and M. Inagaki. Regulatory mechanisms and functions of intermediatefilaments: A study using site- and phosphorylation state-specific antibodies.Cancer Sci., 97:167–174, 2006.13. J. Johnson and K. Omland. Model selection in ecology and evolution. TrendsEcol. Evol., 19:101–108, 2004. 24/294. J.S. Kim, C.-H. Lee, and P. Coulombe. Modeling the self-organization of keratinintermediate filaments. Biophys. J., 99:2748–2755, 2010.15. A. K¨olsch, R. Windoffer, and R. Leube. Actin-dependent dynamics of keratinfilament precursors. Cell Motil. Cytoskel., 66:976–985, 2009.16. A. K¨olsch, R. Windoffer, T. W¨urflinger, T. Aach, and R. Leube. Thekeratin-filament cycle of assembly and disassembly. J. Cell Sci., 123:2266–2272,2010.17. R. Leube, M. Moch, A. Kolsch, and R. Windoffer. “panta rhei”: Perpetualcycling of the keratin cytoskeleton. Bioarchitecture, 1:39–44, 2011.18. Y. Li, A. Brown, and P. Jung. Deciphering the axonal transport kinetics ofneurofilaments using the fluorescence photoactivation pulse-escape method. Phys.Biol., 11:026001, 2014.19. A. Madzvamuse, A.J. Wathen, and P.K. Maini. A moving grid finite elementmethod applied to a model biological pattern generator. J. Comp. Phys.,190:478–500, 2003.20. MATLAB. MATLAB and Statistics Toolbox Release 2013b. The MathWorksInc., Natick, Massachusetts, 2013.21. M. Moch, G. Herberich, T. Aach, R. Leube, and R. Windoffer. Measuring theregulation of keratin filament network dynamics. Proc. Natl. Acad. Sci.,110:10664–10669, 2013.22. C. M¨ohl, N. Kirchgessner, C. Sch¨afer, B. Hoffmann, and R. Merkel. Quantitativemapping of averaged focal adhesion dynamics in migrating cells by shapenormalization. J. Cell Sci., 125:155–165, 2012.23. X. Pan, R. Hobbs, and P. Coulombe. The expanding significance of keratinintermediate filaments in normal and diseased epithelia. Curr. Opin. Cell Biol.,25:47–56, 2013.24. S. Portet and J. Arino. An in vivo intermediate filament assembly model. Math.Biosc. Eng., 6:117–134, 2009.25. S. Portet, O. Arino, J. Vassy, and D. Schoevaert. Organization of the cytokeratinnetwork in an epithelial cell. J. Theor. Biol., 223:313–333, 2003.26. L. Ramms, G. Fabris, R. Windoffer, N. Schwarz, R. Springer, C. Zhou, J. Lazar,S. Stiefel, N. Hersch, U. Schnakenberg, T. Magin, R. Leube, R. Merkel, andB. Hoffmann. Keratins as the main component for the mechanical integrity ofkeratinocytes. Proc. Natl. Acad. Sci., 110:18513–18518, 2013.27. A. Robert, H. Herrmann, M. Davidson, and V. Gelfand. Microtubule-dependenttransport of vimentin filament precursors is regulated by actin and by theconcerted action of rho- and p21-activated kinases. Fased J., 28:2879–2890, 2014.28. K. Seltmann, A. Fritsch, J. Kas, and T. Magin. Keratins significantly contributeto cell stiffness and impact invasive behavior. Proc. Natl. Acad. Sci.,110:18507–18512, 2013.29. N. Snider and B. Omary. Post-translational modifications of intermediatefilament proteins: mechanisms and functions. Nat. Rev., 15:163–177, 2014.25/290. P. Strnad, R. Windoffer, and R.E. Leube. Induction of rapid and reversiblecytokeratin filament network remodeling by inhibition of tyrosine phosphatases.J. Cell Sci., 115:4133–4148, 2002.31. C. Sun, R. Leube, R. Windoffer, and S. Portet. A mathematical model for thekeratin cycle of assembly and disassembly. IMA J. Appl. Math., 80:100–114,2015.32. D. Toivola, P. Strnad, A. Habtezion, and B. Omary. Intermediate filaments takethe heat as stress proteins. Trends Cell Biol., 20:79–91, 2010.33. R. Windoffer, M. Beil, T. Magin, and R. Leube. Cytoskeleton in motion: thedynamics of keratin intermediate filaments in epithelial cells. J. Cell Biol.,194:669–678, 2011.34. R. Windoffer and R. Leube. Detection of cytokeratin dynamics by time-lapsefluorescence microscopy in living cells. J. Cell Sci., 112:4521–4534, 1999.35. S. Woll, R. Windoffer, and R. Leube. Dissection of keratin dynamics: differentcontributions of the actin and microtubule systems. Eur. J. Cell Biol.,84:311–328, 2005.36. K. Yoon, M. Yoon, R. Moir, S. Khuon, F. Flitney, and R. Goldman. Insights intodynamic properties of keratin intermediate filaments in living epithelial cells. J.Cell Biol., 153:503–516, 2001.
Supporting Information Legends
S1 Video
Dynamics of the keratin network in a cell.
Time-lapse fluorescence microscopy ofhepatocellular carcinoma-derived PLC clone PK18-5 stably expressing fluorescent fusionprotein HK18-YFP [30] depicting the dynamic properties of the keratin filaments over atime period of 15 hours. Bar 10 µm . 26/29 − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o1 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o4 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o13 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o16 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o25 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o28 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o2 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o5 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o14 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o17 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o26 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o29 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o3 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o6 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o15 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o18 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o27 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o30 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o7 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o10 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o19 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o22 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o31 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o34 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o8 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o11 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o20 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o23 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o32 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o35 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o9 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o12 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o21 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o24 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o33 − − Lo ca t i on ( µ m ) Concentration ( µ M) . S ce n a r i o36 F i g u r e . B e s t fi t f o r e a c h o f t h e s c e n a r i o s . R e d c u r v e i s t h e r e s p o n s e I i ( x , t f i n a l , ˆ p ) a t h o u r s o f S ce n a r i o i . B l a c k c i r c l e s a r e d a t aa t h o u r s a s i n F i g u r e . . F i g u r e . B e s t f i t f o r e a c h o f t h e s c e n a r i o s . R e d c u r v e i s t h e r e s p o n s e I i ( x , t f i n a l , ˆ p ) a t h o u r s o f S ce n a r i o i . B l a c k c i r c l e s a r e d a t aa t h o u r s a s i n F i g . . .
20 −10 0 10 20010203040 Location ( µ m) C on ce n t r a t i on o f s o l ub l e poo l ( µ M ) −20 −10 0 10 200200400600800 Location ( µ m) C on ce n t r a t i on o f i n s o l ub l e poo l ( µ M )
20 30 40 500.950.9520.9540.9560.9580.960.962 Time (in hours) P r opo r t i on o f I n s o l ub l e −20 −10 0 10 200246810 Location ( µ m) R a t e k ass ( ⋅ )k dis ( ⋅ ) Figure 8. Details about the best model Scenario 21. k ass = 9 . µM/s , k S = 570 . µM , k dis = 0 . µM/s leadingto k max = 1 . µM/s in (16) and k I = 976 . µM/s . Figure 8. Details about the best model Scenario 21. k ass = 9 . µM/s , k S = 570 . µM , k dis = 0 . µM/s leading to k max = 1 . µM/s in (16) and k I = 976 . µM/s . 28/29
20 −10 0 10 2001234 x 10 −3 Location ( µ m) S p ee d ( µ m / s ) Datau(x)
Figure 9. Estimate of the space-dependent speed.
Estimate of thespace-dependent speed given in (13) for the variable speed u ( x ) of the active transportof the insoluble pool. −20 −10 0 10 20Location ( µ m) Data: Sourcek ass (x)
Figure 10. Profile for localized assembly rate of type Sources.
The profile for k ass ( x ) defined in (14) with k max = 1, obtained by fitting the profile of assemblyregions Sources published in [21] and shown in Fig. 2.3. −20 −10 0 10 20Location ( µ m) Data: Sinkk dis (x)
Figure 11. Profile for localized disassembly rate of type Sinks.
The profile for k dis ( x ) defined in (15) and k maxmax