\documentclass{aims}
\usepackage{lineno,hyperref}
\usepackage{caption}
\usepackage{amsmath}
\modulolinenumbers[5]

\usepackage{xcolor}
\usepackage{longtable}
\usepackage{soul}
\usepackage{graphicx}



\newcommand{\stress}{\mbox{\boldmath{$\sigma$}}}
\newcommand{\bD}{\mathbf D}
\newcommand{\diver}{\nabla\cdot}
\newcommand{\Tr}{\mbox{Tr}}
\newcommand{\vpoly}{v_{poly}}
\newcommand{\vperiph}{v_{periph}}
\newcommand{\epsmem}{\varepsilon_{mem}}
\newcommand{\dotepsmem}{\dot{\varepsilon}_{mem}}
\newcommand{\bvarepsilon}{\boldsymbol{\varepsilon}}

\def\typeofarticle{Research article}
\def\currentvolume{20}
\def\currentissue{2}
\def\currentyear{2023}
\def\currentmonth{}
\def\ppages{2408--2438.}
\def\DOI{10.3934/mbe.2023113}
\def\Received{06 August 2022}
\def\Revised{27 October 2022}
\def\Accepted{30 October 2022}
\def\Published{21 November 2022}
\setcounter{page}{1}

\usepackage{cite}
\setcounter{page}{1}

\setlength{\parindent}{2em}
\tolerance=1
\emergencystretch=\maxdimen
\hyphenpenalty=10000
\hbadness=10000
\makeatletter
\renewcommand{\@biblabel}[1]{#1\hfill \hspace{-0.2cm}}
\makeatother

\newcommand{\ep}{\varepsilon}
\newcommand{\eps}[1]{{#1}_{\varepsilon}}



\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\arabic{table}}  


%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\title{On the significance of membrane unfolding in mechanosensitive cell spreading: Its individual and synergystic effects}


\author{% 
	Magdalena A. Stolarska\affil{1,}\corrauth\hspace{0.2em}
	and 
	Aravind R. Rammohan\affil{2}}
	% \shortauthors is used in copyright information in the end of the paper
\shortauthors{the Author(s)}

\address{%
	\addr{\affilnum{1}}{Department of Mathematics, 2115 Summit Ave., University of St. Thomas, St. Paul, MN 55105,
		USA}
	\addr{\affilnum{2}}{Corning Life Sciences, Corning Inc., 836 North St, Tewksbury, MA 01876, USA}}

% corresponding author
\lineskip+\parskip\corraddr{Email: mastolarska@stthomas.edu.}

\maketitle	

\begin{center}
	\textbf{Supplementary}
\end{center}


\section{Justification for additive decomposition of rate of deformation and strain}

The fundamental variable in large strain analysis is the deformation gradient $\mathbf F$, which transforms a material element $d\mathbf X$ in the undeformed or original configuration to the corresponding deformed element $d\mathbf x$ in the deformed or current configuration (see Figure S1A): 
\begin{equation}
d\mathbf x = \mathbf F \, d \mathbf X.
\end{equation}

To incorporate finite active deformations, such as unfolding or exocytosis of the membrane, the multiplicative decomposition of $\mathbf F$ into a passive and active part, $\mathbf F^P$ and $\mathbf F^A$, respectively, is typically used \cite{Lee1969EPD}.  As illustrated in Figure S1A, $\mathbf F^A$ maps $d\mathbf X$ from the original configuration to $d \mathbf x'$ in an intermediate configuration that is assumed to be locally stress-free, whereas $\mathbf F^P$ then maps $d \mathbf x'$ into  $d\mathbf x$ in the current configuration.  The passive deformation results in residual stresses that arise from keeping the deformable body continuous.

In the axisymmetric model of cell spreading described here, deformations of the cell membrane can be described in terms of the radius of the cell (see Figure S1B).  Let $R_0$, $r^A$, and $r$ be the radius of the cell in the original, intermediate, and current configurations, respectively.  In this context, we can define the passive and active components of the deformation gradient in the cell membrane as follows
\begin{eqnarray}
\mathbf F \; = \; \mathbf F^P \mathbf F^A & = & \Big(\underbrace{\frac{r}{r^A}(\mathbf e_{\theta}\otimes \mathbf e_{\theta})}_{\mathbf F^P}\Big)\cdot\Big(\underbrace{\frac{r^A}{r}(\mathbf e_{\theta}\otimes \mathbf e_{\theta})}_{\mathbf F^A}\Big)  \\
& = & \frac{r}{r^A}\frac{r^A}{R_0}(\mathbf e_{\theta}\otimes \mathbf e_{\theta}) \\
& = & \frac{r}{R_0} (\mathbf e_{\theta}\otimes \mathbf e_{\theta}),
\end{eqnarray}
where $\mathbf e_{\theta}$ is the unit vector tangent to the circular membrane.    From this point forward, we drop the dyadic product $(\mathbf e_{\theta}\otimes \mathbf e_{\theta})$ and define the  scalar passive, active, and total components of the deformation gradient in the direction tangent to the membrane as
\begin{equation}
F^P = \frac{r}{r^A} \;\;\;\;\;\;\;\; F^A = \frac{r^A}{R_0}\;\;\;\;\;\;\; F=\frac{r}{R_0}
\label{defgradscalar}
\end{equation}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=4.0in]{FigureS1.png}
\end{center}
\caption{ Decomposition into active and passive components: (A) Conceptual illustration of multiplicative decomposition of the deformation gradient, $\mathbf F$.  $\mathbf F = \mathbf F^P\mathbf F^A$, where the active component ($\mathbf F^A$) maps $d\mathbf X$ from the original into an intermediate, locally stress-free configuration, and the passive component ($\mathbf F^P$) maps the resulting vector to the the current configuration, resulting in $d\mathbf x$.  (B) An illustration of original, intermediate, and passive configurations in the axisymmetric problem considered.  $R$, $r^A$, and $r$, are the cell radii in the original, intermediate, and current configuration, respectively.  $\mathbf e_{\theta}$ is a unit vector tangent to the cell periphery in all three configurations.}
\label{FigS1}
\end{figure}



The Almansi-Hamel strain tensor used in this work is defined as
\begin{equation}
\bvarepsilon \; = \; \frac{1}{2} \Big(1-\mathbf F^{-T}\mathbf F^{-1}\Big).
\end{equation}

In terms of the scalar deformation gradient components defined in (\ref{defgradscalar}), the total, passive, and active strains in the membrane are, respectively, given by
\begin{eqnarray}
\varepsilon_{mem}& = & \frac{1}{2}\Big(1-\Big(\frac{R_0}{r}\Big)^2 \Big) \label{memstrain}\\
\varepsilon^p_{mem} & = & \frac{1}{2}\Big(1-\Big(\frac{r^A}{r}\Big)^2 \Big) \\
\varepsilon^A_{mem} & = & \frac{1}{2}\Big(1-\Big(\frac{R_0}{r^A}\Big)^2 \Big)  \label{memstrainA}.
\end{eqnarray}

Decomposing the inverse of the total deformation in (\ref{memstrain}) into the inverse of the active and passive parts results in
\begin{eqnarray}
\varepsilon_{mem} &  = & \frac{1}{2}\Big(1-\Big(\frac{R_0}{r}\Big)^2 \Big) \; = \; \frac{1}{2}\Bigg( 1 - \Big(\Big(\frac{R_0}{r^A}\Big)\Big(\frac{r^A}{r}\Big)\Big)^2 \Bigg) \\
& = &  \frac{1}{2}\Bigg( 1 - \Big(\frac{r^A}{r}\Big)^2 + \Big(\frac{r^A}{r}\Big)^2 -  \Big(\Big(\frac{R_0}{r^A}\Big)\Big(\frac{r^A}{r}\Big)\Big)^2 \Bigg) \\
& = &   \frac{1}{2}\Bigg( 1 - \Big(\frac{r^A}{r}\Big)^2\Bigg) \; + \; \Big(\frac{r^A}{r}\Big)^2 \cdot \frac{1}{2}\Bigg( 1 - \Big(\frac{R_0}{r^A}\Big)^2\Bigg) \\
& = & \varepsilon^p_{mem} + \Big(\frac{r^A}{r}\Big)^2 \varepsilon^A_{mem} \label{memstrainfinal}
\end{eqnarray}

With the simplifying assumption that passive deformations are small, $r \approx r^A$, which gives $\frac{r^A}{r} \approx 1$.  This approximation reduces (\ref{memstrainfinal}) to an additive decomposition of strains given by
\begin{equation}
\varepsilon_{mem} \; = \;  \varepsilon^p_{mem}  + \varepsilon^A_{mem}.
\end{equation}

Because all active deformations are locally stress-free, stresses in the membrane result entirely from passive deformations.  The passive strain and passive strain rate, given the justification above, are 
\begin{equation}
\varepsilon^p_{mem} = \varepsilon_{mem} - \varepsilon^A_{mem} \;\;\;\;\;\;\;  \dot{\varepsilon}^p_{mem} = \dot{\varepsilon}_{mem} - D^A_{mem}.
\end{equation}

It needs to be noted, however, that the right hand sides of the above equations are fully nonlinear.
$D^A_{mem}$ is the strain rate representing membrane unfolding or exocytosis and needs to be prescribed to determine response of the system.  The algorithmic definition of $D^A_{mem}$ is provided in the main text of this work.  

$D^A_{mem}$ and the polymerization speed, $v_{poly}$, are computed for a given time step, and the velocity of the intracellular actin, $v_c$, is then determined using the computational approach described below.  Via numerical integration this allows us to compute $\varepsilon_{mem}$ and $\varepsilon_{mem}^A$ consistent with their nonlinear definitions specified herein.

The total strain, $\varepsilon_{mem}$ can be determined by integrating its rate consistent with (\ref{memstrain})
\begin{equation}
\dotepsmem \; = \; \frac{R_0^2}{r^3}\big(\vpoly+v_c\big)
\label{memstrainrate}
\end{equation}
(which reproduces Eq (2.9) in the main text.)  To find $\varepsilon^A_{mem}$ we start with the time derivative of $F^A$ as written in (\ref{defgradscalar}):
\begin{equation}
\dot{F}^A \; = \; \frac{\dot{r}^A}{R_0} \; = \; \frac{\dot{r}^A}{r^A} \frac{r^A}{R_0} \; = \; D^A_{mem} \frac{r^A}{R_0} \label{epsmem1} 
\end{equation}
which simplifies to
\begin{equation}
\dot{r^A} \; = \; D^A_{mem}r^A \label{epsmem2}
\end{equation}
and where $D^A_{mem} = \frac{\dot{r}^A}{r^A}$ by definition.  Considering that at a given time $t$ the values of $r^A$ and $D^A_{mem}$ are known (computed), numerical integration of $\dot{r}^A$ gives $r^A(t+\Delta t)$, and $\varepsilon^A_{mem}(t+\Delta t)$ is found using (\ref{memstrainA}).


\section{Parameters used}

The following tables list all parameter values used in the various modeling components.


\begin{table}[H]
	\caption{Parameters used for MM Model.}\setlength{\tabcolsep}{5.9mm}
	\begin{center}
		\label{tableparamsMM}
		\begin{tabular}{llll}
			\hline
			%\multicolumn{4}{|c|} {\em MM Model} \\ \hline
			{Parameters} & {Description} & {Value} & {Reference} \\
			\hline
$K$ & actin bulk viscosity & 16.7 Pa-min & \cite{Kim2014DFB} \\ 
$E_m = E_1+E_2$ & overall membrane elasticity & 10,000 pN & \cite{Hochmuth1973MEM,Vernerey2014MMC} \\ 
$ E_1 $ & membrane elasticity constant & 5000 pN & assumed\\ 
$ E_2 $ & membrane elasticity constant & 5000 pN & assumed \\ 
$\eta_m$ & membrane viscosity & 100 pN-min & \cite{Evans1976MV} \\ 
$\vpoly$ & constant actin polymerization speed & 3 $\mu$m/min &\cite{Mogilner1996CMD,Carlsson2003GVB} \\ 
$\eta_{att}$ & constant attachment drag coefficient & 2 Pa-min/$\mu$m & \cite{Shao2012CAF,Larripa2006T1V} \\ 
  $\alpha_a$ & substrate velocity scale factor & 0.85 & assumed \\ 
  $E_s$ & substrate Young's modulus& 2.5--100 kPa& various \\ 
  $\nu_p$ & substrate Poisson ratio & 0.4 & \cite{Takigawa1996PRP} \\ \hline
  \end{tabular}
  \end{center}
  \end{table}\vspace{-0.7cm}  
  
  \begin{center}
  	\begin{table}[H]
  		\label{tableparamsFA}
  		\caption{Parameters used for FA Model.}\setlength{\tabcolsep}{3.6mm}\vspace{0.25cm}
  		\begin{tabular}{llll}
  			\hline
  			% \multicolumn{4}{|c|} {\em FA Model} \\ \hline
  			{Parameters} & {Description} & {Value} & {Reference} \\
  			\hline
  $k_{on}^0$ & sase FA complex activation rate & 0.0007 (min-$\mu m^2)^{-1}$& assumed \\ 
   $N_{max}$ & total amount of free FA complex & 100 $\mu m^2$ & assumed \\ 
   $\phi_{max}$ & local maximum FA complex concentration & 2.5 a.u. & assumed \\ 
   $k_{off}^0$ & sase FA complex degradation rate & 1.0 $(\mbox{min})^{-1}$ & assumed \\ 
   $k_a$ & attachment traction feedback parameter & 0.4 (Pa)$^{-1}$  & assumed \\ 
   $D$ & bound integrin diffusion coefficient & 0.008 $\mu m^2$/min  & \cite{Lepzelter2010CDI} \\ 
   $\eta_a^0$ & minimum attachment drag coefficient & 0.1 Pa-min/$\mu$m & assumed \\ 
   $\eta_{max}$ & maximum attachment drag & 20 Pa-min/$\mu$m  & \cite{Barnhart2011ADS} \\ 
   $\phi_0$ & attachment force feedback parameter &  0.35 a.u.  & assumed \\ \hline
   \end{tabular}
   \end{table}
    
  \begin{table}[H]
  	\caption{Parameters used for Poly Model.}\setlength{\tabcolsep}{6.9mm}\vspace{0.25cm}
  	\label{tableparamsPoly}
  	\begin{tabular}{llll}
  		\hline
  		% \multicolumn{4}{|c|} {\em Poly Model} \\ \hline
  		{Parameters} & {Description} & {Value} & {Reference} \\
  		\hline
   $\vpoly^0$ & minimum polymerization speed & 0.02 $\mu$m/min & \cite{Welf2013BCB} \\ 
   $\vpoly^{max}$ & maximum polymerization speed & 10.2  $\mu$m/min & \cite{Mogilner1996CMD,Carlsson2003GVB} \\ 
   $ \alpha_v$ & FA complex sensitivity parameter & 0.01 & assumed \\ \hline
   \end{tabular}
   \end{table}
   
  \begin{table}[H]
  	\label{tableparamsExo}
  	\caption{Parameters used for Unfold Model.}\setlength{\tabcolsep}{4.6mm}\vspace{0.25cm}
  	\begin{tabular}{llll}
  		\hline
  		%\multicolumn{4}{|c|} {\em Exo Model} \\ \hline
  		{Parameters} & {Description} & {Value} & {Reference} \\
  		\hline
    $g_{mem}$ & early stage growth rate & 0.01 pN/min & assumed \\ 
    $k_{\phi}$ & FA complex - membrane feedback parameter & 4.0  & assumed \\ 
    $\phi^{TOT}_{max}$ & maximum total FA complex & 100 $\mu m^2$ & assumed \\ 
    $d_m $ & Membrane reserves decay rate & 0.32 (min)$^{-1}$ & assumed \\ \hline
     \end{tabular}
   \end{table}
   
  \begin{table}[H]
 	\label{tableparamsContract}
 	\caption{Parameters used in Contract Model.}\setlength{\tabcolsep}{7.1mm}\vspace{0.25cm}
 	\begin{tabular}{llll}
 		\hline
 		% \multicolumn{4}{|c|} {\em Contract Model} \\ \hline
 		{Parameters} & {Description} & {Value} & {Reference} \\
 		\hline  
    $D^{contract}$ & actomyosin contraction rate & -0.1 (min)$^{-1}$ & \cite{Reymann2012ANA} \\ 
    $k_{\sigma}$ & cellular tension feedback paramter & 10 (Pa)$^{-1}$ & assumed \\ \hline
    
\end{tabular}
\end{table}
\end{center}\vspace{-0.7cm}

The lipid bilayer of the cell membrane is very stiff with small viscosity.    We choose a cell membrane viscosity that is 2--3 orders of magnitude larger than that listed in the literature for the lipid bilayer \cite{Evans1976MV} because such a choice reduces oscillations that may otherwise occur in the numerical method.  We also note that even with a larger membrane viscosity, the elastic stiffness of the membrane dominates the overall material response.  In addition, the actin cortex attached to the cell membrane is more viscous and significantly less stiff than the lipid bilayer and also contributes to the material response at the cell periphery \cite{Berret2016LVL}.

In the FA model, the value of FA complex $\phi(r,t)$ is given in dimensionless arbitrary units (a.u.).  We note that in this particular model component we are not aiming for a quantitative comparison to experimental results, so the parameters for the FA model were chosen so that total FA complex sensitivity to the substrate stiffness qualitatively matches experimental data.   Because $\phi^{TOT}$ is the integral of $\phi(r,t)$ over the cell area, its units are $\mu m^2$.

In the Poly model, the parameters for the minimum and maximum polymerization speeds were chosen from a combination of experimental and modeling results.  We choose the fixed polymerization velocity in the MM model, $\vpoly$, by finding the average value of (2.14) of the main text.  Specifically, we find that given the noted parameter values for $\vpoly^0$ and $\vpoly^{max}$, 
\begin{equation}
\frac{1}{N_{max}}\int_0^{N_{max}}\vpoly(\phi^{TOT})\,d\phi^{TOT} = 3.17
\end{equation}
and therefore chose a value of $\vpoly = 3.0\; \mu m/min$ for the fixed polymerization speed.  

For an arbitrary line segment of length $L$, the active rate of deformation measures $\frac{\Delta L}{L\Delta t}$, where $\Delta L$ is the change in length per time interval $\Delta t$.   Reymann et al. \cite{Reymann2012ANA} measure the speeds of myosin contraction in various actin structures.  From these speeds and the size of actin rings used in their experiments we are able to calculate the contraction rate, $D^{contract}$.  


\section{Computational approach}

We solve the model equations in Matlab R2021a using the finite element method.  Below is an outline of the numerical algorithm that was used to solve the governing equations.  We outline the algorithm for the Contract model since this is the model that includes all mechanisms we consider.  
\begin{enumerate}
\item[1)] {Prior to iterations in time:}
\begin{enumerate}
\item Substrate mesh is created.  For the simulations shown in this work the substrate domain is circular with $radius_{sub} \; \le 50\, \mu$m.  400 radial finite elements are used in the discretization.
\item Cell mesh is created.  For the simulations shown in this work the initial cell radius is $R_0 \le 10\, \mu$m.  100 finite elements are used for the cell throughout the simulation, which employs a moving mesh.  The boundary node initially at $R_0 = 10\,\mu$m  moves with the velocity of the cell periphery, $v_{periph}$.  The remaining nodes are redistributed at each time step so that they remain equally spaced. 
\item {Comment:}  We have also run simulations of the Contract model using cell meshes with 50 and 200 cell elements and have found that the difference in cell spread area sensitivity to substrate stiffness obtained with 100 elements and that obtained with 200 elements is negligible.
\item The initial values of all fields for which we are solving are zero.
\end{enumerate}
\item[2)] {Within each iteration in time:} Simulations correspond to a time interval of 30 min.  For all simulation results shown here we use a time step $\Delta t = 0.005$ min.
\begin{enumerate}
\item Rate of active deformation representing actomyosin contractions, $D^{contract}$, as used in (2.23) of the main text, is defined on the current mesh.  $D^{contract}$ is set to the value specified in Table S5.  This constant rate of contraction acts from the cell center to a distance 2 $\mu$m from the cell periphery.  For all other models in which actomyosin contraction is not present, $D^{contract}$ is set to zero.
\item Given $v_c(r,t)$, $\sigma_r(r,t)$, $\sigma_\theta(r,t)$, $u_s(r,t)$, $\phi(r,t)$, $\phi^{TOT}(t)$, and $v_s(r,t)$ (i.e. the substrate velocity), we use the finite element method to solve for $v_c(r,t+\Delta t)$, $\sigma_r(r,t+\Delta t)$, $\sigma_\theta(r,t+\Delta t)$ and $\Delta u_s(x,t+\Delta t)$ (i.e. the change is substrate displacement over time $\Delta t$).
\begin{enumerate}
\item The polymerization velocity, $v_{poly}$ is calculated using (2.14) in the main text.
\item The finite element method requires that the governing equations are put into a weak form and discretized in time.  We use an implicit time stepping method to calculate updated values of $v_c$, $\sigma_r$, $\sigma_\theta$ and $\Delta u_s$. The time-discretized weak forms of governing equations are given below.  For clarity we have omitted explicit indication of which variables depend on $r$.
$$
\mbox{}
$$
{Weak form of (2.3) from main text:}
\begin{multline}
\int_{-1}^1 \Big(\frac{\partial \overline{v}}{\partial r}  \sigma_r(t+\Delta t) \, + \, \overline{v} \frac{1}{r}\sigma_\theta(t + \Delta t)\Big)rJ\,d\xi \, + \\
\hspace{-0.5cm}\int_{-1}^1 \overline{v}\eta_{\text{att}}(\phi(t)) \big(v_c(t+\Delta t) - v_s(t)\big) rJ\, d\xi \,  = \big(\overline{v}\sigma_{mem}(t)\big)\Big|_{r_N}
\label{weakeqcell}
\end{multline}
In this and all weak forms below, the domain of integration is a reference element, which is parameterized by $-1\le \xi \le 1$, where the spatial coordinate local to the reference element is $\xi$.  The variable $J$ represents the Jacobian of transformation between global radial coordinates in the physical space, denoted by $r$, to the space of local coordinates, $\xi$.  For all integrals we use Gaussian quadrature with two quadrature points to integrate over a reference element.    In (\ref{weakeqcell}), the variable $\overline{v}$, is the test function.  The attachment viscosity, $\eta_{\text{att}}$ is calculated in terms of the local value of $\phi(r,t)$ according to (2.13) in the main text.  Lastly, $r_N$ denotes the position of the boundary node at the cell periphery on which the membrane stress, $\sigma_{mem}(t)$ is acting.  Calculation of $\sigma_{mem}$ is described further below.
$$
\mbox{}
$$
{Weak form of (2.4) and (2.5) from main text:}
\begin{eqnarray}
\hspace{-1cm}\int_{-1}^1 \overline{\sigma}_r\Big( \sigma_r(t+\Delta t) - K\big(\frac{\partial v_c}{\partial r}(t+\Delta t) - D^{contract} \big) \Big) rJ\,d\xi & = & 0 \label{weaksr}\\
\hspace{-1cm}\int_{-1}^1 \overline{\sigma}_\theta\Big(\sigma_\theta(t+\Delta t) - K \frac{v_c(t+\Delta t)}{r}\Big) r J \,d\xi & = & 0 \label{weakst}
\end{eqnarray}
In both (\ref{weaksr}) and (\ref{weakst}), we have explicitly set $\lambda = 0$, since this is the value used in simulations.  The variables $\overline{\sigma}_r$ and $\overline{\sigma}_\theta$ are test functions.\\
\\
{Weak form of (2.11) from main text:}
\begin{multline}
\int_{-1}^1 \Bigg[C_s\frac{\partial\overline{u}_s}{\partial r}\frac{\partial{\Delta u_s}}{\partial r}(t+\Delta t) \, +\\
 \hspace{-1cm} C_s\nu_p \Big(\frac{\partial{\overline{u}_s}}{\partial r}\frac{\Delta u_s(t+\Delta t)}{r} +  \frac{\overline{u}_s}{r}\frac{\partial \Delta u_s}{\partial r}(t+\Delta t)\Big)  +  C_s \frac{\overline{u}_s}{r}\frac{\Delta u_s(t+\Delta t)}{r}\, +\\
   \eta_{\text{att}}(\phi(t))\Big( v_s(t) - 
 v_c(t+\Delta t) \Big)  \Bigg]rJ\, d\xi  \, = \\
\hspace{-1cm}- \int_{-1}^1 \Bigg[C_s\frac{\partial\overline{u}_s}{\partial r}\frac{\partial{u_s}}{\partial r}(t) \, + 
 C_s\nu_p \Big(\frac{\partial{\overline{u}_s}}{\partial r}\frac{u_s(t))}{r} +  \frac{\overline{u}_s}{r}\frac{\partial u_s}{\partial r}(t)\Big)  +  C_s \frac{\overline{u}_s}{r}\frac{u_s(t)}{r}\Bigg]r J\, d\xi
 \label{weaksub}
\end{multline}
In (\ref{weaksub}), $\overline{u}_s$ is a test function.
\item The weak forms in (\ref{weakeqcell})--(\ref{weaksub}) are spatially discretized.  We use $C^0$ piecewise-linear Lagrange finite elements to discretize the cell velocity, $v_c$, the change in substrate displacement $\Delta u_s$, and the test function $\overline{v}$ and $\overline{u}$.  We use piecewise constant finite elements to discretize the radial and angular cell stresses, $\sigma_r$ and $\sigma_\theta$, respectively, and test functions  $\overline{\sigma}_r$ and $\overline{\sigma}_\theta$.     In order to prevent spurious oscillations is stresses, it is typical to use finite element spaces of lower order to discretize stresses in mixed finite element methods \cite{Noor1983SAS}.  The fully discretized governing equations are assembled into a global matrix with the following structure
\begin{equation}
\hspace{-2.1cm}\left[\begin{array}{cccc} M_{\overline{v}v_c} & M_{\overline{v}\sigma_r} & M_{\overline{v}\sigma_\theta}  & 0 \\
M_{\overline{\sigma}_rv_c} & M_{\overline{\sigma}_r\sigma_r} & 0 & 0  \\ 
M_{\overline{\sigma}_\theta v_c}  & 0 & M_{\overline{\sigma}_\theta\sigma_\theta} & 0 \\
M_{\overline{u}v_c}& 0 & 0 & M_{\overline{u}\Delta u_s} 
\end{array}\right]
\left[\begin{array}{c} v_c(t+\Delta t) \\ \sigma_r(t+\Delta t) \\ \sigma_\theta(t+\Delta t) \\ \Delta u_s(t+\Delta t) \end{array}\right] \; = \; \left[\begin{array}{c} R_{v_c} \\ R_{\sigma_r} \\ R_{\sigma_t} \\ R_{\Delta u_s} \end{array}\right] 
\label{arraysolve}
\end{equation}
In this linear system, the variables denoted by $M$ are assembled matrices corresponding to components of the weak forms above.  Similarly, variables $R$ correspond to the components of the weak form evaluated at time $t$ or those to those that are constant.  Using (\ref{arraysolve}), we solve for $v_c(x,t+\Delta t)$, $\sigma_r(x,t+\Delta t)$, $\sigma_{\theta}(t+\Delta t)$ and $\Delta u_s(x,t+\Delta t)$ at all nodal positions.
\end{enumerate}
\item Given the current value of $v_c(x,t+\Delta t)$ at the boundary node corresponding to the cell periphery, the velocity of the cell periphery, $v_{periph}$, is updated, and the position of the boundary node at the cell periphery, $r_N$, is updated according to $r_N (t+\Delta t) = r_N(t) + v_{periph}\Delta t$.  The positions of the remaining nodes are updated so that they are equally spaced.  We record the mesh velocity, which is given by
\begin{equation}
V_{mesh}(r_i) = \frac{v_{periph}}{N_{el}}i
\label{vmesh}
\end{equation}
where $r_i$ is the position of the $i^{th}$ node and $N_{el}$ is the number of elements.
\item Given the current value of $v_c(x,t+\Delta t)$ at the boundary node corresponding to the cell periphery, the updated mesh, and active rate of membrane deformation at time $t$, we update the active rate of membrane deformation $D^A_{mem}(t+\Delta t)$ and membrane tension $\sigma_{mem}(t+\Delta t)$. 
\begin{enumerate}
\item To update the membrane tension and begin the algorithm (for $t\ge0.7$ min.) for calculating the rate of active membrane deformation as described in Section 2.5 of the main text, we compute the updated rate of change of the membrane strain, $\dot{\varepsilon}_{mem}(t+\Delta t)$ according to (\ref{memstrainrate}).  Similarly, the total membrane strain, $\varepsilon_{mem}(t+\Delta t)$,  and the total active strain, $\varepsilon^A_{mem}(t+\Delta t)$ are updated according to (\ref{memstrain}) and (\ref{memstrainA}).  Given these values, one can calculate the change in membrane stress, $\Delta \sigma_{mem}(t+\Delta t)$ from the discretized version of (2.17) of the main text.  This discretized version is given by
\begin{multline}
\frac{\Delta \sigma_{mem}(t+\Delta t)}{\Delta t} = (E_1 + E_2)\big(\dot{\varepsilon}_{mem}(t+\Delta t) - D^A_{mem}(t)\big) + \\
\hspace{-1cm}\frac{E_1E_2}{\eta_m}\big(\varepsilon_{mem}(t+\Delta t) - \varepsilon^A_{mem}(t+\Delta t)\big) - \frac{E_2}{\eta_m}\big(\sigma_{mem}(t) + \Delta\sigma_{mem}(t+\Delta t)\big)
\label{deltasig}
\end{multline}
We calculate the updated membrane tension from 
\begin{equation}
\sigma_{mem}(t+\Delta t) \, = \, \sigma_{mem}(t) + \Delta \sigma_{mem}(t+\Delta t)
\end{equation}
\item The next step of the algorithm in Section 2.5 of the main text requires reducing the rate of membrane tension in a manner proportional to $\phi^{TOT}(t)$, as given by equation (2.19) of the main text.  The discrete form of this equation is given by
\begin{equation}
\Delta\sigma_{mem}^{D^A} = \, e^{-k_{\phi}f(\phi^{TOT}(t))}\Delta{\sigma}_{mem}(t+\Delta t)
\end{equation}
where $\Delta{\sigma}_{mem}(t+\Delta t)$ is determined from (\ref{deltasig}).
\item The reduced change in membrane tension is substituted into (\ref{deltasig}) to compute the rate of active membrane deformation that would result in this value of change in membrane tension.  Specifically, the following equation is solved for $\overline{D}^A_{mem}$.
\begin{multline}
\frac{\Delta \sigma_{mem}^{D^A}}{\Delta t} = (E_1 + E_2)\big(\dot{\varepsilon}_{mem}(t+\Delta t) - \overline{D}^A_{mem}\big) + \\
\hspace{-1cm}\frac{E_1E_2}{\eta_m}\big(\varepsilon_{mem}(t+\Delta t) - \varepsilon^A_{mem}(t+\Delta t)\big) - \frac{E_2}{\eta_m}\sigma_{mem}(t+\Delta t)
\label{deltasigMEM}
\end{multline}
\item Once  $\overline{D}^A_{mem}$ is calculated, this value is reduced in time, as given by (2.21) in the main text, to obtain an updated active rate of membrane deformation to be used at the next time step.
\end{enumerate}
\item Given $\phi^{TOT}(t)$, $\sigma_r(t+\Delta t)$, $v_c(t+\Delta t)$, and a substrate velocity updated according to (2.7) of the main text, $v_s(t+\Delta t)$, we use the finite element method to update the spatial distribution of the FA complex, $\phi(r,t+\Delta t)$.
\begin{enumerate}
\item The attachment force $f_{\text{attach}}$ is updated according to (2.6) of the main text.
\item { Comment}: We use an Arbitrary Eulerian-Lagrangian (ALE) approach because the finite element mesh for the cell is updated at each time step. Given a general mesh velocity, $\mathbf V_{mesh}$, a general ALE formulation for the reaction-diffusion-advection of scalar quantity $c$ is given by the following (see  \cite{MacDonald2016CMC}):
\begin{equation}
\frac{\partial c}{\partial t} + \nabla\cdot\big(c\mathbf v\big) - \mathbf V_{mesh}\nabla c \, = \, D\nabla^2 c +R(c)
\end{equation}
where $\mathbf v$ is the advection velocity, $D$ is the diffusion coefficient, and $R(c)$ describes any reaction terms.  A general overview of ALE approaches can be found in \cite{Donea2004ALE}.
\item  The numerical solution to our model equations requires that the evolution equation for $\phi$ accounts for the mesh velocity.  In our case, the mesh velocity is denoted by $V_{mesh}$ and computed via (\ref{vmesh}).  Note that since $\phi$ represents the entire activated FA complex and is localized within the cell via attachment to the substrate, we do not assume that $\phi$ is advected with the velocity of the actin network.  The temporally discretized weak form of (2.24) of the main text, along with the addition of the contribution from the moving mesh, is given by
\begin{multline}
\int_{-1}^1 \overline{\phi}\Bigg(\frac{\phi(t+\Delta t) - \phi(t)}{\Delta t} - V_{mesh}\frac{\partial\phi}{\partial r}(t+\Delta t) + \\ 
\hspace{-1.5cm}k_{off}^0 e^{-k_a|f_{attach}(t+\Delta t)|}e^{-k_\sigma h(\sigma_r(t+\Delta t))}\phi(t+\Delta t)\Bigg)rJ\,d\xi  \, + \, \int_{-1}^1 D\frac{\partial\overline{\phi}}{\partial r} \frac{\partial \phi}{\partial r}(t+\Delta t) rJ\, d\xi \, = \\
\int_{-1}^1\overline{\phi} k_{on}^0\big(N_{max}-\phi^{TOT}(t)\big)\big(\phi_{max} - \phi(t)\big)rJ\,d\xi
\end{multline}
Again, explicit writing of $r$ dependence has been omitted for clarity.  Here, $\overline{\phi}$ are test functions.  We use $C^0$ piecewise linear finite elements for the spatial discretization of $\phi$ and $\overline{\phi}$.
\item Lastly, $\phi^{TOT}(t+\Delta t)$ is computed by integrating $\phi(r,t+\Delta t)$ over the current cell area:
\begin{equation}
\phi^{TOT}(t+\Delta t) \, = \, 2\pi \int_0^{r_N} \phi(r,t+\Delta t) r \,dr
\end{equation}
\end{enumerate}
\end{enumerate}
\end{enumerate}



\label{computation}
\pagebreak[4]
\begin{figure}[H]
	\begin{center}
		\includegraphics[width=5in]{FigureS2.png}
	\end{center}\vspace{0.2cm}
	\caption{ Effect of changing parameter $\lambda$.
		In our simulations we assume that $\lambda=0$. (See Eqs (2.3) and (2.4) of the main text). Changing the parameter $\lambda$ does not significantly affect quantitative results and has no affect on the qualitative conclusions of this work.  To quantify the affect of changing $\lambda$ from $\lambda=0$ to $\lambda=5$, we calculate the percent increase in cell spread  in all models we consider in this work.  For the MM model, the percent increase in cell spread area for $\lambda = 0$ is 2.50\% and for $\lambda = 5$ it is 2.56\%.   For the FA model, the percent increase in spread area for $\lambda = 0$ is 16.66\% and for $\lambda = 5$ it is 17.00\%.  For the Poly model, the percent increase in spread area for $\lambda = 0$ is 38.40\%. and for $\lambda = 5$ it is 39.49\%.  For the Unfold model, the percent increase in spread area for $\lambda = 0$ is 267.62\% and for $\lambda = 5$ it is 290.89\%.   For the Contract model, the percent increase in spread area for $\lambda = 0$ is 203.96\% and for $\lambda = 5$ it is 247.19\%.  }
\end{figure}
\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=5in]{FigureS3.png}
	\end{center}\vspace{0.2cm}
	\caption{Substrate velocities. Substrate velocities computed for a substrate Young's modulus $E_s = 100$ kPa for each of the five models we consider.  Due to actin retrograde flow, substrate displacements, and therefore substrate velocities, are oriented toward cell interior.  Maximum substrate speeds occur near the cell periphery. }
\end{figure}
\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=5.5in]{FigureS4.png}
	\end{center}\vspace{0.2cm}
	\caption{A comparison of attachment tractions.
		In all figures white arrows indicate direction of traction that the cell exerts onto the substrate.  (A) Attachment tractions for MM model.  (B) Attachment tractions for FA model.  (C) Experimentally measured attachment tractions from \cite{Oakes2014GRT} (reproduced with permission).  Cells were plated on prepatterned polyacrylamide gels and imaged after 6--24 hours of spreading.  Fibronectin patterns controlled for cell spread area, which in these images is 800 $\mu m^2$.  Color bar for experimental attachment tractions is found in left-most panel. }
\end{figure}

\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=5.in]{FigureS5.png}%
	\end{center}\vspace{0.2cm}
	\caption{Comparison of polymerization speeds, retrograde flow speeds at the cell periphery, velocities of the cell periphery, and membrane tension for different models.  Row 1 shows actin polymerization speeds.  Row 2 shows the retrograde flow speeds.  Row 3 shows the velocity of the cell periphery, which is the (positive) polymerization velocity added to the (negative) retrograde flow velocity.  Row 4 is the membrane tension.  The figures in column 1 correspond to data for the Const Poly - Unfold model.  Column 2 is for the Poly model.  Column 3 is for the Unfold model. }
\end{figure}
\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=4.6in]{FigureS6.png}
	\end{center}\vspace{0.15cm}
	\caption{Membrane unfolding / exocytosis active deformation. (A) The total membrane area added by unfolding / exocytosis. 
		(B)
		Time evolution of the active strain rate, $D_{mem}^A$, representing membrane  unfolding / exocytosis.    The inset shows the time evolution of the membrane unfolding rate for times $t = 1$ to 4.5 min, where largest differences in active rate of deformation occur. }
\end{figure}
\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=4.6in]{FigureS7.png}
	\end{center}\vspace{0.15cm}
	\caption{A comparison of the effects of membrane unfolding and adhesion dependent polymerization speed with contractility on equilibrium cell spread area: Top left image shows spread area dependence on substrate Young’s modulus for the Contant Poly / Unfold model  with contractility.   In this case, the equilibrium cell spread area increases by 60.96\% as the substrate Young’s modulus increases from 2.5 kPa to 100 kPa.  The bottom left image shows spread area dependence on substrate Young’s modulus for the Poly model with contractility.  Here, the equilibrium cell spread area increases by 37.88\% as the substrate Young’s modulus increases from 2.5 kPa to 100 kPa.  When both FA dependent actin polymerization speed and membrane unfolding / exocytosis are combined with contractility (i.e., the Contraction model)  the two mechanisms are amplified and the equilibrium cell area increases by 203.96\%.  In general, the combination of membrane unfolding / exocytosis with FA-dependent actin polymerization amplifies the dependence of cell spread area on substrate stiffness.  The addition of contractility slightly decreases the percent increase in spread area for all combinations of mechanisms.}
\end{figure}


\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=4.5in]{FigureS8.png}
	\end{center}\vspace{0.15cm}
	\caption{Difference in retrograde flow speed of Contract model and Unfold model. Quantity graphed is (retrograde flow speed for Contract model) - (retrograde flow speed for Unfold model).  In (A) a substrate stiffness of $E_s = 2.5$ kPa is used, and in (B) a substrate stiffness of $E_s = 100$ kPa.  Vertical dashed line indicates transition region at which retrograde flow speed of Contract model exceeds that of Unfold model.  Arrow indicates region where retrograde flow speed of Contract model is larger.  Horizontal axis represents finite element (FEM) node id.  We chose this as the horizontal variable instead of spatial coordinates, which change as the cell spreads, so that we can compare differences in retrograde flow speeds in the same position {\em relative to} the cell center.}
\end{figure}



\newpage
\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=4.6in]{FigureS9.png}
	\end{center}\vspace{0.15cm}
	\caption{Changing value of $k_{\sigma}$ in Contract model. (A) Cell spread area dependence on substrate stiffness for Contract model with $k_{\sigma} = 100$ (Pa)$^{-1}$.  The spread area increases by 193.44\%, which is slightly smaller compared to results for $k_{\sigma} = 10$ (Pa)$^{-1}$shown in Figure~7 of main text. (B)  Intracellular radial stresses at $t=30$ min for Contract model with $k_{\sigma} = 100$ (Pa)$^{-1}$ and substrate stiffness $E_s = 2.5\, kPa$ (top) and $E_s = 100\, kPa$ (bottom). Intracellular radial stresses are almost identical to those for Contract model with $k_{\sigma} = 10$ (Pa)$^{-1}$ shown in Figure 7 of main text. Black circles indicate regions of intracellular tension. (C) Comparison of sensitivity of total FA complex to substrate stiffness for the FA model, Unfold model, and Contract model. All three models (and others not illustrated) follow the same trend. }
\end{figure}

\begin{thebibliography}{999} 
  
\bibitem{Lee1969EPD}
E.~H. Lee, Elastic-plastic deformation at finite strains, {\it{ASME J. Appl. Mech.}}, {\bf 36} (1969), 1--6.
\doilink{https://doi.org/10.1115/1.3564580}

\bibitem{Kim2014DFB}
T.~Kim, M.~L. Gardel, E.~Munro, Determinants of fluidlike behavior and
effective viscosity in cross-linked actin networks, {\it{Biophys. J.}},
{\bf 106} (2014), 526--534.
\doilink{https://doi.org/10.1016/j.bpj.2013.12.031}

\bibitem{Hochmuth1973MEM}
R.~M. Hochmuth, N.~Mohandas, P.~Blackshear~Jr, Measurement of the elastic
modulus for red cell membrane using a fluid mechanical technique, {\it{Biophys.
		J.}}, {\bf 13} (1973), 747--762.
\doilink{https://doi.org/10.1016/S0006-3495(73)86021-7}

\bibitem{Vernerey2014MMC}
F.~J. Vernerey, M.~Farsad, A mathematical model of the coupled mechanisms of
cell adhesion, contraction and spreading, {\it{J. Math. Biol.}},
{\bf 68} (2014), 989--1022.
\doilink{https://doi.org/10.1007/s00285-013-0656-8}

\bibitem{Evans1976MV}
E.~Evans, R.~Hochmuth, Membrane viscoelasticity, {\it{ Biophys. J.}}, {\bf 16}
(1976), 1--11.
\doilink{https://doi.org/10.1016/S0006-3495(76)85658-5}

\bibitem{Mogilner1996CMD}
A.~Mogilner, G.~Oster, Cell motility driven by actin polymerization,
{\it{Biophys. J.}}, {\bf 71} (1996), 3030--3045.
\doilink{https://doi.org/10.1016/S0006-3495(96)79496-1}

\bibitem{Carlsson2003GVB}
A.~Carlsson, Growth velocities of branched actin networks, {\it{Biophys. J.}},
{\bf 84} (2003), 2907--2918.
\doilink{https://doi.org/10.1016/S0006-3495(03)70018-6}

\bibitem{Shao2012CAF}
D.~Shao, H.~Levine, W. J. Rappel, Coupling actin flow, adhesion, and morphology
in a computational cell motility model, {\it{PNAS}}, {\bf 109} (2012), 6851--6856.
\doilink{https://doi.org/10.1073/pnas.1203252109}

\bibitem{Larripa2006T1V}
K.~Larripa, A.~Mogilner, Transport of a 1d viscoelastic actin--myosin strip of
gel as a model of a crawling cell, {\it{Physica A}}, {\bf 372} (2006), 113--123.
\doilink{https://doi.org/10.1016/j.physa.2006.05.008}

\bibitem{Takigawa1996PRP}
T.~Takigawa, Y.~Morino, K.~Urayama, T.~Masuda, Poisson's ratio of
polyacrylamide (paam) gels, {\it{Polym. Gels Networks}}, {\bf 4} (1996), 1--5.
\doilink{https://doi.org/10.1016/0966-7822(95)00013-5}

\bibitem{Lepzelter2010CDI}
D.~Lepzelter, M.~H. Zaman, Clustered diffusion of integrins, {\it{Biophys.
		J.}}, {\bf 99} (2010), L106--L108.
\doilink{https://doi.org/10.1016/j.bpj.2010.11.007}

\bibitem{Barnhart2011ADS}
E.~L. Barnhart, K. C. Lee, K.~Keren, A.~Mogilner, J.~A. Theriot, An
adhesion-dependent switch between mechanisms that determine motile cell
shape, {\it{PLoS Biol.}}, {\bf 9} (2011), e1001059.
\doilink{https://doi.org/10.1371/journal.pbio.1001059}

\bibitem{Welf2013BCB}
E.~S. Welf, H.~E. Johnson, J.~M. Haugh, Bidirectional coupling between
integrin-mediated signaling and actomyosin mechanics explains
matrix-dependent intermittency of leading-edge motility, {\it{Mol. Biol.
		Cell}}, {\bf 24} (2013), 3945--3955.
\doilink{https://doi.org/10.1091/mbc.E13-06-0311}

\bibitem{Reymann2012ANA}
A. C. Reymann, R.~Boujemaa-Paterski, J. L. Martiel, C.~Gu{\'e}rin, W.~Cao,
H.~F. Chin, et al., Actin network
architecture can determine myosin motor activity, {\it{Science}}, {\bf 336} (2012),
1310--1314.
\doilink{https://doi.org/10.1126/science.1221708}

\bibitem{Berret2016LVL}
J.~F. Berret, {Local viscoelasticity of living cells measured by rotational
	magnetic spectroscopy}, {\it{Nat. Commun.}}, {\bf 7} (2016), 1--9.
\doilink{https://doi.org/10.1038/ncomms10134}

\bibitem{Noor1983SAS}
A. K.~Noor, W. D. Pilkey, \emph{State-of-the-Art Surveys on Finite Element Technology},  American Society of Mechanical Engineers,  New York, 1983.

\bibitem{MacDonald2016CMC}
G.~MacDonald, J.~A. Mackenzie, M.~Nolan, R.~Insall, A computational method for
the coupled solution of reaction--diffusion equations on evolving domains and
manifolds: Application to a model of cell migration and chemotaxis, {\it{J. Comput. Phys.}}, {\bf 309} (2016), 207--226.
\doilink{https://doi.org/10.1016/j.jcp.2015.12.038}

\bibitem{Donea2004ALE}
J.~Donea, A.~Huerta, J. P. Ponthot, A.~Rodríguez-Ferran, {\it{Arbitrary
		Lagrangian–Eulerian Methods}}, John Wiley \& Sons, Ltd, (2004), 413--437.
\doilink{https://doi.org/10.1002/0470091355.ecm009}

\bibitem{Oakes2014GRT}
P.~W. Oakes, S.~Banerjee, M.~C. Marchetti, M.~L. Gardel, Geometry regulates
traction stresses in adherent cells, {\it{Biophys. J.}}, {\bf 107} (2014),
825--833.
\doilink{https://doi.org/10.1016/j.bpj.2014.06.045}
\end{thebibliography}


\end{document}