September 25, 2017, Monday

# Introduction

A critical question in biology is to understand the mechanisms by which biological tissues are able to adaptively change their properties and structure in response to external stimuli. Indeed these mechanisms are at the origin of important phenomena such as morphognesis, remodeling and growth, wound healing and cancer dynamics. Various studies have shown that tissue adaptation was driven by the ability of contractile cells to apply forces and deform their surrounding extra-cellular matrix. Cells acquire these abilities through the polymerization of actin and myosin monomers into a network of contractile fibers called Stress Fibers (SF). These SF connect the cells to their surrounding Extracellular Matrix (ECM) via the focal adhesions (FA,multi-molecular complexes localized on the cells' membrane) and allow them to sense stimuli conveyed by their ECM and respond accordingly by changing their internal stiffness and contractility. Studies have shown that SF organization was strongly influenced by the nature of the cells' mechanical environment (such as stiffness, deformation, loading frequency). When tissue are not stretch, it has been shown that cell contractility and stress fiber density was strongly dependent on substrate stiffness; i.e. the stiffer the subtrate the stronger the contraction. Furthermore, in the presence of stretch, SF were shown to align in the direction of maximum stretch for static loading and perpendicular with the direction of stretch for cyclic loading.

Here we propose a mathematical model based on thermodynamics to predict orientation and contractile stress in SF in response to uniaxial constant or cyclic stretch or changes in the substrate properties (such as stiffness, Poisson's ratio, anisotropy). Using the mechano-chemical potential of monomers and SF polymers to find the angular volume fraction of SF when the system reaches chemical equilibrium, the model shows the effects of substrate stiffness on the SF organization, and can capture the cell's behavior under both constant and cyclic stretch with great consistency with experimental observations. This model is first solved numerically, than analytically and we then finally address the double question of under what loading context the model is most sensitive to what parameter. This will give us an idea of how to grow contractile cells with the highest volume fraction of SF.

(click on the pictures to see them normally)

# The Model

The capacity of contractile cells like fibroblasts to exert tensional forces and change their internal stiffness is principally due to the existence of a network of SF that generates tension through acto-myosin cross bridges dynamics. }. In vitro experiments have shown that SF are generated by self-assembly when the cell is subjected to a stiff substrate or to constant or cyclic stretch, each case showing different SF organization patterns In order to capture those different SF organizations, a model is constructed in which cells are deposited on a substrate of stiffness $E_{s}$ with maximum confluence. Each cell maintain a pool of monomers (G-actin) that can be recruited to be polymerized into SF (F-actin). This polymerization reaction is viewed as a thermodynamic phenomenon dependent on the difference of mechano-chemical potential between monomers and polymers and on the angle from the direction of stretch. The final distribution of SF is obtained when both mechanical and chemical equilibria are reached for all $\theta \in \left[0,\pi \right]$.

Mechanical equilibrium. A cell at mechanical equilibrium in direction $\theta$ is subjected to three stresses: a stress $\sigma _{{\theta }}^{s}$ from the deformation of the substrate to which it is attached, a stress $\sigma _{{\theta }}^{p}$ from the passive deformation of its own and neighbor cells' cytoskeleton and a third stress from its internal contraction $\sigma _{{\theta }}^{c}$ (Fig.2.). This yields the following equilibrium:

$\sigma _{{\theta }}^{p}+\sigma _{{\theta }}^{s}+\sigma _{{\theta }}^{c}=0$(1)

Chemical equilibrium. The polymerization of G-actin into F-actin driving the self-assembly of stress fibers is described using the chemical potential of monomers $\mu ^{m}$ and polymers $\mu ^{p}$. The chemical equilibrium is obtained when minimizing the Gibbs free energy \cite{PeterOUP1985}:

${\frac {\partial G}{\partial \xi }}=\sum (\nu _{i}\mu _{i})=\mu ^{p}-\mu ^{m}=0$(2)

With $\xi$ the extent of the reaction, and $\nu _{i}$ and $\mu _{i}$ respectively the stoichiometric coefficients and chemical potentials (here $\mu ^{p}$ and $\mu ^{m}$) of the different constituents involved in the reaction. If $\mu ^{p}>\mu ^{m}$, stress fibers are disassembled into monomers, if $\mu ^{p}<\mu ^{m}$ monomers are polymerized into stress fibers, and the state of equilibrium is reached when the potential of monomers and polymers are equal.

Constitutive equations. Here, we derive a simple mathematical model that can capture many features of SF organization and contraction based on the following three ideas:

• Contractile stress stabilizes and promotes the stress fibers self-assembly and the elastic energy of a stress fiber causes its disassembly. Indeed, it has been shown that together with the RhoA enzyme, the contractility of the cell promotes the formation of SF \cite{ChrzanowskaRUP1996}: the activation of RhoA drives the myosin filament assembly and its interaction with actin, and the contractile tension allows the alignment and consolidation of the different components \cite{walcottPNAS2010}. Furthermore, myofribrils, which are close in structure to SF, start fragmenting as their strain , i.e. their elastic energy increases \cite{SimpsonCR1999}. The same behavior has been observed for SF \cite{LeeBBRC2010}. We can now construct the chemical potentials $\mu ^{m}$ and $\mu _{{\theta }}^{p}$ of SF building unit. It consists in the classical formulation of the chemical potential plus the mechano-chemical potential of the SF \cite{PathakJTRSI2008,ShemeshPNAS2005}:

$\mu ^{m}=\mu _{{0}}^{m}+kT_{{temp}}\ln \left({\frac {c^{m}}{c_{{0}}}}\right)$(3)

$\mu _{{\theta }}^{p}=\mu _{{0}}^{p}+kT_{{temp}}\ln \left({\frac {c_{{\theta }}^{p}}{c_{{0}}}}\right)+V_{{CU}}\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle$(4)

where $c^{m}$ and $c_{{p}}^{{\theta }}$ are the concentrations of monomers and polymers, $\mu _{{0}}^{{m,p}}$ and $c_{{0}}$ are the reference chemical potentials and concentration of monomers and polymers, k is the Boltzmann constant, $T_{{temp}}$ the absolute temperature and$V_{{CU}}$ is the volume of SF contractile unit. The average volumic mechanical contribution $\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle$ to the polymer's chemical potential on one stretching period $T$ is constructed as follows:

$\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle ={\frac {1}{T}}\int _{{0}}^{{T}}\left({\frac {1}{2}}E(t)\varepsilon _{{\theta }}^{{2}}-\sigma ^{{*}}\varepsilon _{{\theta }}\right)dt$(5)

The elastic energy term ${\frac {1}{2}}E(t)\varepsilon _{{\theta }}^{{2}}$ will increase the polymer's chemical potential, contributing to the disassembly of stress fibers, and the contractile work term $\sigma ^{{*}}\varepsilon _{{\theta }}$ will decrease it, inducing a stabilization of stress fibers.

• Stress fibers exhibit a visco-elastic behavior. It has been reported that SF behave like a viscoelastic cable,tensed by the action of actomyosin motors \cite{kumarBJ2006} (fig1.c). The time-dependent elasticity of a SF is therefore constructed as follows:

$E(t)=E_{1}+E_{2}\exp(-t/\tau )$(6)

With $E_{1}$ and $E_{2}\exp(-t/\tau )$ respectively the SF linear elasticity and time dependent elasticity, $t$ the relaxation time and $\tau$ the relaxation time constant for a SF.

• The cell contractily and stiffness increase proportionally with SF density in given directions. Further assuming that both cell and substrate behave in a linear elastic fashion, the mechanical constitutive equations and the resulting equilibrium equation as follows:

$\sigma _{{\theta }}^{s}=E_{s}\varepsilon _{{\theta }}$(7)

$\sigma _{{\theta }}^{p}=(E_{c}+\Phi _{{\theta }}^{p}E(t))\varepsilon _{{\theta }}$(8)

$\sigma _{{\theta }}^{c}=\Phi _{{\theta }}^{p}\sigma ^{{*}}$(9)

where $E_{s}$ and $E_{c}$ are respectively the stiffnesses of substrate cell's cytoskeleton and stress fibers, $\varepsilon _{{\theta }}$ is the strain in direction $\theta$, $\sigma ^{{*}}$ is the stress produced by one stress fiber and $\Phi _{{\theta }}^{p}$ is the volume fraction of stress fibers in direction $\theta$.

# Solution

First to be solved is the mechanical equilibrium for either stress free or displacement control boundary conditions to get the angular strain $\varepsilon _{{\theta }}$. This allows us to calculate the average mechanical contribution to the polymer's potential $\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle$ and subsequently solve the chemical equilibrium:

$\mu _{{0}}^{m}+kT_{{temp}}\ln \left({\frac {c^{m}}{c_{{0}}^{m}}}\right)=\mu _{{0}}^{p}+kT_{{temp}}\ln \left({\frac {c_{{\theta }}^{p}}{c_{{0}}^{p}}}\right)+V_{{CU}}\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle$(10)

Using the volume fraction-concentration relation ships $\Phi _{{\theta }}^{{p}}=(M^{a}/\rho ^{a})c_{{\theta }}^{{p}}$ and $\Phi ^{{m}}=(M^{a}/\rho ^{a})c^{{m}}$ and the conservation of mass of polymers and monomers $\Phi ^{m}+\left\langle \Phi _{{\theta }}^{{p}}\right\rangle =1$, one can solve equation (10) for the angular volume fraction of polymers:

$\Phi _{{\theta }}^{{p}}=\Phi ^{m}\exp {\left(-{\frac {1}{KT_{{temp}}}}\left(\mu _{{0}}^{{p}}-\mu _{{0}}^{{m}}+V_{{CU}}\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle \right)\right)}$(11)

$\Phi ^{m}=1-{\frac {1}{\pi }}\int _{{0}}^{{\pi }}\Phi _{{\theta }}^{{p}}d\theta$(12)

# Numerical solving

We first solve the set of equation above using Matlab's nonlinear equation solver on the discretized region $\left[0,\pi \right]$:

• effect of substrate stiffness:

Figure a shows the effect of substrate stiffness on the level of production (volume fraction) of SF. Figure b shows the angular volume fraction for a isotropic and anisotropic substrate.

• effect of constant stretch:

Figure a shows the effect of stretch on the volume fraction of SF. We can see that passed a critical stretch, the volume fraction of SF starts decreasing. Figure b shows the angular volume fraction for different level of stretch.

• effect of cyclic stretch:

Figures a and b show the effect of cyclic stretch for Poisson's ratio of respectively 0 and 0.5.

For the effect of substrate stiffness, the system can only be solved numerically, but from the numerical results, it appears clear that the volume fraction of SF remains relatively low, and that the volume fraction of monomers stays almost constant. We can therefore try to find an analytical solution to this problem with the approximation of low SF volume fraction and constant monomers volume fraction. As for the constant and cyclic stretch cases, the solution are already analytical.

# Analytical solution for the effect of substrate stiffness case

Here is the equation for the effect of substrate's stiffness case: $\Phi _{{\theta }}^{{p}}=\Phi ^{{m}}\exp \left(-{\frac {1}{KT}}\left(\mu _{{0}}^{{p}}-\mu _{{0}}^{{m}}+V_{{cu}}\left(\underbrace {{\frac {1}{2}}E_{{1}}({\frac {\Phi _{{\theta }}^{{p}}\sigma ^{{*}}}{E_{{p}}}})^{2}}_{{{\mbox{negligible}}}}+{\frac {(\Phi _{{\theta }}^{{p}})(\sigma ^{{*}})^{2}}{E_{{p}}}}\right)\right)\right)$(13)

With $E_{1}$ and $E_{{p}}$ being of the same order, we can now neglect the term in $(\Phi _{{\theta }}^{{p}})^{2}$ with the approximation of SF volume fraction remaining low, and monomer volume fraction $\Phi ^{m}$ constant. We now have an equation of the form:

$x=a\exp(-b-c*x)$(14)

solvable using the product log function, or Lambert function W of order 0:

$\Phi _{{\theta }}^{{p}}={\frac {W\left(\Phi ^{m}e^{{-1/KT(\mu _{{0}}^{{p}}-\mu _{{0}}^{{m}})}}{\frac {V_{{cu}}(\sigma ^{{*}})^{2}}{KTE_{{p}}}}\right)}{{\frac {V_{{cu}}(\sigma ^{{*}})^{2}}{KTE_{{p}}}}}}$(15)

Equation (15) is the analytical solution for the angular volume fraction of SF in the effect of substrate's stiffness case. Here is a comparison between analytical (green) and numerical (blue)solution:

They seem to match quite well with a slight offset in the high stiffness region.

The analytical solutions for the constant and cyclic stretch cases are given by:

$\Phi _{{\theta }}^{{p}}=\Phi ^{m}\exp {\left(-{\frac {1}{KT_{{temp}}}}\left(\mu _{{0}}^{{p}}-\mu _{{0}}^{{m}}+V_{{CU}}\left\langle \mu _{{mech}}^{{p}}(\varepsilon _{{\theta }})\right\rangle \right)\right)}$(11)

with $\Phi ^{m}$ constant and

$\varepsilon _{{\theta }}=\varepsilon _{{11}}\left((1+\nu )\cos ^{2}{(\theta )}-\nu \right)$(16)

with $\varepsilon _{{11}}$ a given input in the constant stretch case, and in the cyclic stretch case

$\varepsilon _{{11}}=\varepsilon _{{0}}^{{c}}+\varepsilon ^{{c}}\cos \left({\frac {2\pi }{T}}t\right)\;\;{\mbox{with}}\;\;\varepsilon _{{0}}^{{c}}=0\;\;{\mbox{and}}\;\;\varepsilon ^{{c}}=0.1$(17)

# Sensitivity Analysis

With the analytical solutions for the three cases (effect of substrate's stiffness, effect of constant and cyclic stretching) we can now perform a sensitivity analysis. We chose the Elementary Effect Method. we want to investigate the effect of the biological parameters of the cell such as the contractile stress of SF, the different stiffnesses of SF and cytoskeleton, as well as the experiment parameters such as the substrate Poisson's ratio and the stretching frequency in the three different contexts of (1) effect of substrate's stiffness , (2) effect of constant stretching, (3) and effect of cyclic stretching to find what what parameter in what context can the most influence the growth of SF.

To do so, we analyse the effect $d_{i}$ of a small change successively on each paramater $p_{i}$ on the function $\Phi ^{p}$:

$d_{i}({\mathbf {p}}^{j})={\frac {\Phi ^{p}(p_{{1}},p_{2},...,p_{{i-1}},p_{{i}}+\Delta ,p_{{i+1}},...,p_{k})}{\Delta }}$

with the different ${\mathbf {p}}^{j}$ being r different random sets of sampling points obtained with a Latin Hypercube Square design.

From these elementary effects, we can then calculate the mean and variance for each parameters, for the 3 different contexts:

mean: $\mu _{i}={\frac {1}{r}}\sum _{{j=1}}^{{r}}d_{i}({\mathbf {p}}^{j})$

variance: $\sigma _{i}={\sqrt {{\frac {1}{r-1}}\sum _{{j=1}}^{{r}}(d_{i}({\mathbf {p}}^{j})-\mu _{i})^{2}}}$

The mean represents the overall influence of the parameter on the volume fraction, and the variance estimates the impact of the parameter due to interaction with the others. Both these values need to be high for a parameter to have a real impact. Here we plot the mean and variance for each parameters in the three cases considered:

• case 1: effect of substrate: sensitivity of the volume fraction of SF to the contractile stress of SF$\sigma ^{{*}}$ and cytoskeleton stiffness $E_{c}$:

• case 2: effect of constant stretch: we add the sensitivity to the Poisson's ratio (since now it is not isotropic anymore) and the SF's siffness $E_{1}$:

• case 3: effect of cyclic stretch: with the sensitivity to the time dependent part of the SF's stiffness $E_{2}$ and the frequency of stretching $\delta$:

# Conclusion

From this sensitivity analysis, we can see that the growth of contractile cells is the most affected by the SF level of contractile stress $\sigma ^{{*}}$, and that this sensitivity is higher when the cell is cyclically stretched. It is also worth noticing a rather unexpected result being that the growth of SF is more sensitive to the cytoskelton's stiffness in the constant stretch case than in the cyclic stretch case. These results where obtained for parameters around their commonly accepted values.

# Matlab code

clear all

KT=4.2e-21; mu=17.1e-21; V=1e-23; sigma=10; Ec=1000; epsilon0c=0.1; epsilonc=0.1; theta=0; nu=0.25; E1=100; E2=100; delta=1;

k=2; p=4; d=2/3; r=10,

for kk=1:3

   if kk==1


X=lhsdesign(r,2); X(:,1)=X(:,1)*sigma*2; X(:,2)=X(:,2)*Ec*2;

   for i=1:length(X)
par1(1)=KT;
par1(2)=mu;
par1(3)=V;
par1(4)=X(i,1)+d*X(i,1);
par1(5)=X(i,2);
par1(6)=epsilon0c;
par1(7)=epsilonc;
par1(8)=theta;
par1(9)=nu;
par1(10)=E1;
par1(11)=E2;
par1(12)=delta;

par2(1)=KT;
par2(2)=mu;
par2(3)=V;
par2(4)=X(i,1);
par2(5)=X(i,2)+d*X(i,2);
par2(6)=epsilon0c;
par2(7)=epsilonc;
par2(8)=theta;
par2(9)=nu;
par2(10)=E1;
par2(11)=E2;
par2(12)=delta;

par0(1)=KT;
par0(2)=mu;
par0(3)=V;
par0(4)=X(i,1);
par0(5)=X(i,2);
par0(6)=epsilon0c;
par0(7)=epsilonc;
par0(8)=theta;
par0(9)=nu;
par0(10)=E1;
par0(11)=E2;
par0(12)=delta;

[yS0 yC0 yCy0] = W( 1,par0);
[yS1 yC1 yCy1] = W( 1,par1);
[yS2 yC2 yCy2] = W( 1,par2);
d1S(i)=(yS1-yS0)/d;
d2S(i)=(yS2-yS0)/d;
end



mu1S=1/r*sum(d1S); mu2S=1/r*sum(d2S);

sigma1S=sqrt((1/(r-1))*sum((d1S-ones(1,length(d1S))*mu1S).^2)); sigma2S=sqrt((1/(r-1))*sum((d2S-ones(1,length(d2S))*mu2S).^2));

muS=[mu1S mu2S]; sigmaS=[sigma1S sigma2S]; figure plot(muS',sigmaS','o'); xlabel('mean'); ylabel('variance');

text(muS(1),sigmaS(1),'sigma*')
text(muS(2),sigmaS(2),'Ec')

   elseif kk==2

for i=1:length(X)


X=lhsdesign(r,4); X(:,1)=X(:,1)*sigma*2; X(:,2)=X(:,2)*Ec*2; X(:,3)=X(:,3)*nu*2; X(:,4)=X(:,4)*E1*2;

   par1(1)=KT;
par1(2)=mu;
par1(3)=V;
par1(4)=X(i,1)+d*X(i,1);
par1(5)=X(i,2);
par1(6)=epsilon0c;
par1(7)=epsilonc;
par1(8)=theta;
par1(9)=X(i,3);
par1(10)=X(i,4);
par1(11)=E2;
par1(12)=delta;

par2(1)=KT;
par2(2)=mu;
par2(3)=V;
par2(4)=X(i,1);
par2(5)=X(i,2)+d*X(i,2);
par2(6)=epsilon0c;
par2(7)=epsilonc;
par2(8)=theta;
par2(9)=X(i,3);
par2(10)=X(i,4);
par2(11)=E2;
par2(12)=delta;

par3(1)=KT;
par3(2)=mu;
par3(3)=V;
par3(4)=X(i,1);
par3(5)=X(i,2);
par3(6)=epsilon0c;
par3(7)=epsilonc;
par3(8)=theta;
par3(9)=X(i,3)+d*X(i,3);
par3(10)=X(i,4);
par3(11)=E2;
par3(12)=delta;

par4(1)=KT;
par4(2)=mu;
par4(3)=V;
par4(4)=X(i,1);
par4(5)=X(i,2);
par4(6)=epsilon0c;
par4(7)=epsilonc;
par4(8)=theta;
par4(9)=X(i,3);
par4(10)=X(i,4)+d*X(i,4);
par4(11)=E2;
par4(12)=delta;

   par0(1)=KT;
par0(2)=mu;
par0(3)=V;
par0(4)=X(i,1);
par0(5)=X(i,2);
par0(6)=epsilon0c;
par0(7)=epsilonc;
par0(8)=theta;
par0(9)=X(i,3);
par0(10)=X(i,4);
par0(11)=E2;
par0(12)=delta;

[yS0 yC0 yCy0] = W( 50,par0);
[yS1 yC1 yCy1] = W( 50,par1);
[yS2 yC2 yCy2] = W( 50,par2);
[yS3 yC4 yCy3] = W( 50,par3);
[yS4 yC3 yCy4] = W( 50,par4);
d1C(i)=(yC1-yC0)/d;
d2C(i)=(yC2-yC0)/d;
d3C(i)=(yC3-yC0)/d;
d4C(i)=(yC4-yC0)/d;
end



mu1C=1/r*sum(d1C); mu2C=1/r*sum(d2C); mu3C=1/r*sum(d3C); mu4C=1/r*sum(d4C);

sigma1C=sqrt((1/(r-1))*sum((d1C-ones(1,length(d1C))*mu1C).^2)); sigma2C=sqrt((1/(r-1))*sum((d2C-ones(1,length(d2C))*mu2C).^2)); sigma3C=sqrt((1/(r-1))*sum((d3C-ones(1,length(d3C))*mu3C).^2)); sigma4C=sqrt((1/(r-1))*sum((d4C-ones(1,length(d4C))*mu4C).^2));

muC=[mu1C mu2C mu3C mu4C]; sigmaC=[sigma1C sigma2C sigma3C sigma4C]; figure plot(muC',sigmaC','o'); xlabel('mean'); ylabel('variance');

text(muC(1),sigmaC(1),'sigma*')
text(muC(2),sigmaC(2),'Ec')
text(muC(3),sigmaC(3),'nu')
text(muC(4),sigmaC(4),'E1')


   elseif kk==3

for i=1:length(X)


X=lhsdesign(r,6); X(:,1)=X(:,1)*sigma*2; X(:,2)=X(:,2)*Ec*2; X(:,3)=X(:,3)*nu*2; X(:,4)=X(:,4)*E1*2; X(:,5)=X(:,5)*E2*2; X(:,6)=X(:,5)*delta*2;

   par1(1)=KT;
par1(2)=mu;
par1(3)=V;
par1(4)=X(i,1)+d*X(i,1);
par1(5)=X(i,2);
par1(6)=epsilon0c;
par1(7)=epsilonc;
par1(8)=theta;
par1(9)=X(i,3);
par1(10)=X(i,4);
par1(11)=X(i,5);
par1(12)=X(i,6);

par2(1)=KT;
par2(2)=mu;
par2(3)=V;
par2(4)=X(i,1);
par2(5)=X(i,2)+d*X(i,2);
par2(6)=epsilon0c;
par2(7)=epsilonc;
par2(8)=theta;
par2(9)=X(i,3);
par2(10)=X(i,4);
par2(11)=X(i,5);
par2(12)=X(i,6);

par3(1)=KT;
par3(2)=mu;
par3(3)=V;
par3(4)=X(i,1);
par3(5)=X(i,2);
par3(6)=epsilon0c;
par3(7)=epsilonc;
par3(8)=theta;
par3(9)=X(i,3)+d*X(i,3);
par3(10)=X(i,4);
par3(11)=X(i,5);
par3(12)=X(i,6);

par4(1)=KT;
par4(2)=mu;
par4(3)=V;
par4(4)=X(i,1);
par4(5)=X(i,2);
par4(6)=epsilon0c;
par4(7)=epsilonc;
par4(8)=theta;
par4(9)=X(i,3);
par4(10)=X(i,4)+d*X(i,4);
par4(11)=X(i,5);
par4(12)=X(i,6);

       par5(1)=KT;
par5(2)=mu;
par5(3)=V;
par5(4)=X(i,1);
par5(5)=X(i,2);
par5(6)=epsilon0c;
par5(7)=epsilonc;
par5(8)=theta;
par5(9)=X(i,3);
par5(10)=X(i,4);
par5(11)=X(i,5)+d*X(i,5);
par5(12)=X(i,6);

par6(1)=KT;
par6(2)=mu;
par6(3)=V;
par6(4)=X(i,1);
par6(5)=X(i,2);
par6(6)=epsilon0c;
par6(7)=epsilonc;
par6(8)=theta;
par6(9)=X(i,3);
par6(10)=X(i,4);
par6(11)=X(i,5);
par6(12)=X(i,6)+d*X(i,6);

par0(1)=KT;
par0(2)=mu;
par0(3)=V;
par0(4)=X(i,1);
par0(5)=X(i,2);
par0(6)=epsilon0c;
par0(7)=epsilonc;
par0(8)=theta;
par0(9)=X(i,3);
par0(10)=X(i,4);
par0(11)=X(i,5);
par0(12)=X(i,6);

[yS0 yC0 yCy0] = W( 50,par0);
[yS1 yC1 yCy1] = W( 50,par1);
[yS2 yC2 yCy2] = W( 50,par2);
[yS3 yC4 yCy3] = W( 50,par3);
[yS4 yC3 yCy4] = W( 50,par4);
[yS5 yC5 yCy5] = W( 50,par5);
[yS6 yC6 yCy6] = W( 50,par6);
d1Cy(i)=(yCy1-yCy0)/d;
d2Cy(i)=(yCy2-yCy0)/d;
d3Cy(i)=(yCy3-yCy0)/d;
d4Cy(i)=(yCy4-yCy0)/d;
d5Cy(i)=(yCy5-yCy0)/d;
d6Cy(i)=(yCy6-yCy0)/d;

end



mu1Cy=1/r*sum(d1Cy); mu2Cy=1/r*sum(d2Cy); mu3Cy=1/r*sum(d3Cy); mu4Cy=1/r*sum(d4Cy); mu5Cy=1/r*sum(d6Cy); mu6Cy=1/r*sum(d5Cy);

sigma1Cy=sqrt((1/(r-1))*sum((d1Cy-ones(1,length(d1Cy))*mu1Cy).^2)); sigma2Cy=sqrt((1/(r-1))*sum((d2Cy-ones(1,length(d2Cy))*mu2Cy).^2)); sigma3Cy=sqrt((1/(r-1))*sum((d3Cy-ones(1,length(d3Cy))*mu3Cy).^2)); sigma4Cy=sqrt((1/(r-1))*sum((d4Cy-ones(1,length(d4Cy))*mu4Cy).^2)); sigma5Cy=sqrt((1/(r-1))*sum((d5Cy-ones(1,length(d5Cy))*mu5Cy).^2)); sigma6Cy=sqrt((1/(r-1))*sum((d6Cy-ones(1,length(d6Cy))*mu6Cy).^2));

muCy=[mu1Cy mu2Cy mu3Cy mu4Cy mu5Cy mu6Cy]; sigmaCy=[sigma1Cy sigma2Cy sigma3Cy sigma4Cy sigma5Cy sigma6Cy]; figure plot(muCy',sigmaCy','o'); xlabel('mean'); ylabel('variance');

text(muCy(1),sigmaCy(1),'sigma*')
text(muCy(2),sigmaCy(2),'Ec')
text(muCy(3),sigmaCy(3),'nu')
text(muCy(4),sigmaCy(4),'E1')
text(muCy(5),sigmaCy(5),'E2')
text(muCy(6),sigmaCy(6),'delta')
end


end

function [yS yC yCy] = W( x,par) % KT=4.2e-21; % mu=17.1e-21; % V=1e-23; % sigma=10; % Ec=1000; KT=par(1); mu=par(2); V=par(3); sigma=par(4); Ec=par(5); epsilon0c=par(6); epsilonc=par(7); theta=par(8); nu=par(9); E1=par(10); E2=par(11); delta=par(12);

kappa=(1+nu)*cos(theta)^2-nu;

   C1c=0.5*epsilon0c^2*kappa^2*(E1)-(sigma)*(epsilon0c*kappa);

C1cy=0.5*epsilon0c^2*kappa^2*(E1+E2*exp(-delta)/delta)...
+epsilon0c*epsilonc*kappa^2*(E1+E2*(delta*(1-exp(-delta)))/(4*pi^2+delta^2))...
+0.25*epsilonc^2*kappa^2*(E1+E2*(1-exp(-delta))*(8*pi^2+delta^2)/(16*pi^2*delta+delta^3))...
-(sigma)*(epsilon0c*kappa+epsilonc*kappa);


yS=lambertw(1*exp((-1/(KT))*(mu))*(Ec*sigma^2*V)/(x*KT))/((Ec*sigma^2*V)/(x*KT)); yC=1*exp(-1/KT*(mu+Ec*V*(C1c))); yCy=1*exp(-1/KT*(mu+Ec*V*(C1cy)));

% z=(1+(A*B/x)-sqrt(1+2*A*B/x))/((A^2*B)/(x^2));

% n=0; % for i=1:0.1:100 % n=n+1; % [WW1(n),WW2(n)]=W(i,par); % end % plot([1:0.1:100],WW1,[1:0.1:100],WW2) end