Functions

SetDirectory[NotebookDirectory[]]

bVerbose=True;

numDigitsPrecision=50;

setModelParams[]:=Block[{},(
adjParams=Map[#[[1]]&,adjParamsInit];
unaccountedParams=Complement[allModelParams,Join[Map[#[[1]]&,fixedParams],adjParams]];
nlmParams=allModelParams/.Map[Rule[#[[1]],#[[2]]]&,fixedParams];
If[bVerbose,If[Length[unaccountedParams]==0,Print["All model parameters accounted for."],Print["WARNING! Some model parameters are not accounted for."]]];
If[Length[unaccountedParams]==0,TableForm[N[Join[Map[{#[[1]],"fixed",#[[2]]}&,fixedParams],Map[{#[[1]],"adjustable",#[[2]]}&,adjParamsInit]]],TableHeadings->{None,{"param","fixed/adj","initial value"}}],TableForm[N[Join[Map[{#[[1]],"fixed",#[[2]]}&,fixedParams],Map[{#[[1]],"adjustable",#[[2]]}&,adjParamsInit],Map[{#,"?","?"}&,unaccountedParams]]],TableHeadings->{None,{"param","fixed/adj","initial value"}}]]
)];

Clear[modelERexitLog4,logkϵ,logkm1,logk2,logka,logkb,logkc,logkd,logKE,logKG,k1,km2,kγ,t];
modelERexitLog4[logkϵ_?NumberQ,logkm1_?NumberQ,logk2_?NumberQ,logka_?NumberQ,logkb_?NumberQ,logkc_?NumberQ,logkd_?NumberQ,logKE_?NumberQ,logKG_?NumberQ] :=  (modelERexitLog4[logkϵ,logkm1,logk2,logka,logkb,logkc,logkd,logKE,logKG]=Module[{tmax,AT,A,X,y,ϵ0,ϵ1,ϵ2,γ0,γ1,γ2,p,eqns,initialConditions,system,sol},(
tmax=36;
AT={{-10^logkϵ,0,0,0,0,0,0},{10^logkϵ,-k1,10^logkm1,0,0,0,0},{0,k1,-(10^logkm1+10^logka),0,0,10^logkb,0},{0,0,0,-kγ,0,0,0},{0,0,10^logka,kγ,-(10^logk2+10^logkc),km2,0},{0,0,0,0,10^logk2,-(km2+10^logkb+10^logkd),0},{0,0,0,0,10^logkc,10^logkd,0}};
A=AT/.k1->10^logKE*10^logkm1/.km2->10^logKG*10^logk2/.kγ->0;
X[t_] = {ϵ0[t],ϵ1[t],ϵ2[t],γ0[t],γ1[t],γ2[t],p[t]};
eqns = MapThread[#1==#2&,{X'[t],A.X[t]}];
initialConditions={ϵ0[0]==1,ϵ1[0]==0,ϵ2[0]==0,γ0[0]==0,γ1[0]==0,γ2[0]==0,p[0]==0};
system=Join[eqns,initialConditions];
sol=NDSolve[system,{ϵ0,ϵ1,ϵ2,γ0,γ1,γ2,p},{t,0,tmax}];
Function[{t},First[ϵ0/.sol][t]+First[ϵ1/.sol][t]+First[ϵ2/.sol][t]]
)]);

Clear[modelGolgiExitLog4,logkγ,logkm1,logk2,logka,logkb,logkc,logkd,logKE,logKG,k1,km2,kϵ,t];
modelGolgiExitLog4[logkγ_?NumberQ,logkm1_?NumberQ,logk2_?NumberQ,logka_?NumberQ,logkb_?NumberQ,logkc_?NumberQ,logkd_?NumberQ,logKE_?NumberQ,logKG_?NumberQ] :=  (modelGolgiExitLog4[logkγ,logkm1,logk2,logka,logkb,logkc,logkd,logKE,logKG]=Module[{tmax,AT,A,X,y,ϵ0,ϵ1,ϵ2,γ0,γ1,γ2,p,eqns,initialConditions,system,sol},(
tmax=36;
AT={{-kϵ,0,0,0,0,0,0},{kϵ,-k1,10^logkm1,0,0,0,0},{0,k1,-(10^logkm1+10^logka),0,0,10^logkb,0},{0,0,0,-10^logkγ,0,0,0},{0,0,10^logka,10^logkγ,-(10^logk2+10^logkc),km2,0},{0,0,0,0,10^logk2,-(km2+10^logkb+10^logkd),0},{0,0,0,0,10^logkc,10^logkd,0}};
A=AT/.k1->10^logKE*10^logkm1/.km2->10^logKG*10^logk2/.kϵ->0;
X[t_] = {ϵ0[t],ϵ1[t],ϵ2[t],γ0[t],γ1[t],γ2[t],p[t]};
eqns = MapThread[#1==#2&,{X'[t],A.X[t]}];
initialConditions={ϵ0[0]==0,ϵ1[0]==0,ϵ2[0]==0,γ0[0]==1,γ1[0]==0,γ2[0]==0,p[0]==0};
system=Join[eqns,initialConditions];
sol=NDSolve[system,{ϵ0,ϵ1,ϵ2,γ0,γ1,γ2,p},{t,0,tmax}];
Function[{t},First[γ0/.sol][t]+First[γ1/.sol][t]+First[γ2/.sol][t]]
)]);

allModelParamsLogPS4j={logkϵ,logkγ,logkm1Full,logkm1tmd,logKELAT,logKEAllL,logKEtmd,logka,logkb,logk2,logKGLAT,logKGAllL,logkc,logkd};

modelAllJointLogPS4j[logkϵ_?NumberQ,logkγ_?NumberQ,logkm1Full_?NumberQ,logkm1tmd_?NumberQ,logKELAT_?NumberQ,logKEAllL_?NumberQ,logKEtmd_?NumberQ,logka_?NumberQ,logkb_?NumberQ,logk2_?NumberQ,logKGLAT_?NumberQ,logKGAllL_?NumberQ,logkc_?NumberQ,logkd_?NumberQ,intSpecies_?NumberQ,intCompartment_?NumberQ,time_?NumberQ] :=  (modelAllJointLogPS4j[logkϵ,logkγ,logkm1Full,logkm1tmd,logKELAT,logKEAllL,logKEtmd,logka,logkb,logk2,logKGLAT,logKGAllL,logkc,logkd,intSpecies,intCompartment,time]=Piecewise[{{modelERexitLog4[logkϵ,logkm1Full,logk2,logka,logkb,logkc,logkd,logKELAT,logKGLAT][time],intSpecies==1∧intCompartment==1},{modelGolgiExitLog4[logkγ,logkm1Full,logk2,logka,logkb,logkc,logkd,logKELAT,logKGLAT][time],intSpecies==1∧intCompartment==2},{modelERexitLog4[logkϵ,logkm1tmd,logk2,logka,logkb,logkc,logkd,logKEtmd,logKGLAT][time],intSpecies==2∧intCompartment==1},{modelGolgiExitLog4[logkγ,logkm1tmd,logk2,logka,logkb,logkc,logkd,logKEtmd,logKGLAT][time],intSpecies==2∧intCompartment==2},{modelERexitLog4[logkϵ,logkm1Full,logk2,logka,logkb,logkc,logkd,logKEAllL,logKGAllL][time],intSpecies==3∧intCompartment==1},{modelGolgiExitLog4[logkγ,logkm1Full,logk2,logka,logkb,logkc,logkd,logKEAllL,logKGAllL][time],intSpecies==3∧intCompartment==2},{modelERexitLog4[logkϵ,logkm1tmd,logk2,logka,logkb,logkc,logkd,logKEtmd,logKGAllL][time],intSpecies==4∧intCompartment==1},{modelGolgiExitLog4[logkγ,logkm1tmd,logk2,logka,logkb,logkc,logkd,logKEtmd,logKGAllL][time],intSpecies==4∧intCompartment==2}}]);

plotModelPrediction[]:=Block[{},(
modelCurves=modelCurvesT=dataListCur=dataCur={};
Do[
dataListCur=dataListT[[intSpeciesCur]];
Do[
dataCur=dataListCur[[intCompartmentCur]];
{tmin,tmax}=Map[#[[1]]&,dataCur[[{1,-1}]]];
timeVals=Range[tmin,tmax,(tmax-tmin)/100];
modelCurveCur=Map[{#,Apply[modelCur,Join[N[allModelParams/.Map[Rule[#[[1]],#[[2]]]&,fixedParams]/.Map[Rule[#[[1]],#[[2]]]&,adjParamsInit]],{intSpeciesCur,intCompartmentCur,#}]]}&,timeVals];
modelCurvesT=Append[modelCurvesT,modelCurveCur],
{intCompartmentCur,Length[dataListCur]}];
modelCurves=Append[modelCurves,modelCurvesT];
modelCurvesT={},
{intSpeciesCur,Length[dataListT]}];
legendModel={"ER fit","Golgi fit"};
pmodel=Table[ListPlot[modelCurves[[i]],Joined->True,PlotStyle->color,PlotTheme->{"Scientific","BoldColor","SansLabels"},FrameLabel->{"time [h]","Residence fraction"},PlotLegends->Placed[legendModel,{0.75,0.75}],PlotLabel->plotLabel[[i]],ImageSize->300,AspectRatio->1],{i,Length[modelCurves]}];
Magnify[GraphicsGrid[Partition[MapThread[Show[#1,#2]&,{pdataNoLines,pmodel}],2],ImageSize->{800,800}],0.9]
)];

plotModelFit[]:=Block[{},(
modelCurves=modelCurvesT=dataListCur=dataCur={};
Do[
dataListCur=dataListT[[intSpeciesCur]];
Do[
dataCur=dataListCur[[intCompartmentCur]];
{tmin,tmax}=Map[#[[1]]&,dataCur[[{1,-1}]]];
timeVals=Range[tmin,tmax,(tmax-tmin)/100];
modelCurveCur=Map[{#,Apply[modelCur,Join[N[allModelParams/.Map[Rule[#[[1]],#[[2]]]&,fixedParams]/.bfp],{intSpeciesCur,intCompartmentCur,#}]]}&,timeVals];
modelCurvesT=Append[modelCurvesT,modelCurveCur],
{intCompartmentCur,Length[dataListCur]}];
modelCurves=Append[modelCurves,modelCurvesT];
modelCurvesT={},
{intSpeciesCur,Length[dataListT]}];
pmodel=Table[ListPlot[modelCurves[[i]],Joined->True,PlotStyle->color,PlotTheme->{"Scientific","BoldColor","SansLabels"},FrameLabel->{"time [h]","Residence fraction"},PlotLegends->Placed[legendModel,{0.75,0.75}],PlotLabel->plotLabel[[i]],ImageSize->300,AspectRatio->1],{i,Length[modelCurves]}];
Magnify[GraphicsGrid[Partition[MapThread[Show[#1,#2]&,{pdataNoLines,pmodel}],2],ImageSize->{800,800}],0.9]
)];

example data

(* data are a list of values {time (h), residence fraction} *)
erLAT={{0.,1.},{0.2,0.75},{0.4,0.56},{0.6,0.42},{0.8,0.31},{1.,0.23},{1.2,0.18},{1.4,0.13},{1.6,0.1},{1.8,0.075},{2.,0.056},{2.2,0.042},{2.4,0.032},{2.6,0.024},{2.8,0.018}};

erLATtmd={{0.,1.},{2.,0.80},{4.,0.64},{6.,0.51},{8.,0.41},{10.,0.33},{12.,0.27},{14.,0.21},{16.,0.17},{18.,0.13},{20.,0.11},{22.,0.09},{24.,0.07},{26.,0.057},{28.,0.046},{30.,0.037},{32.,0.03},{34.,0.024}};

erAllL={{0.,1.},{0.60,0.78},{1.2,0.61},{1.8,0.47},{2.4,0.37},{3.,0.29},{3.6,0.23},{4.2,0.17},{4.8,0.14},{5.4,0.11},{6.,0.08},{6.6,0.065},{7.2,0.051},{7.8,0.04},{8.4,0.03},{9.,0.024},{9.6,0.02}};

erAllLtmd={{0.,1.},{2.,0.83},{4.,0.683`},{6.,0.56},{8.,0.47},{10.,0.39},{12.,0.32},{14.,0.26},{16.,0.22},{18.,0.18},{20.,0.15},{22.,0.12},{24.,0.10},{26.,0.084},{28.,0.07},{30.,0.06},{32.,0.04},{34.,0.039}};

ListPlot[{erLAT,erAllL,erLATtmd,erAllLtmd},Joined->True,PlotRange->All,PlotTheme->{"Scientific","BoldColor","SansLabels","OpenMarkers"},FrameLabel->{"time [h]","Residence fraction"},PlotLegends->{"LAT","All-L","LAT-TMD","All-L-TMD"}]

golgiLAT={{0.,1.},{0.4,0.79},{0.8,0.63},{1.2,0.49},{1.6,0.39},{2.,0.31},{2.4,0.25},{2.8,0.19},{3.2,0.15},{3.6,0.12},{4.,0.096},{4.4,0.076},{4.8,0.06},{5.2,0.047},{5.6,0.037}};

golgiLATtmd={{0.,1.},{0.4,0.78},{0.8,0.61},{1.2,0.47},{1.6,0.37},{2.,0.29},{2.4,0.23},{2.8,0.17},{3.2,0.14},{3.6,0.11},{4.,0.083},{4.4,0.065},{4.8,0.05},{5.2,0.04},{5.6,0.03}};

golgiAllL={{0.,1.},{0.8,0.79},{1.6,0.63},{2.4,0.50},{3.2,0.40},{4.,0.32},{4.8,0.25},{5.6,0.20},{6.4,0.16},{7.2,0.13},{8.,0.10},{8.8,0.081},{9.6,0.064},{10.4,0.051},{11.2,0.04}};

golgiAllLtmd={{0.,1.},{0.8,0.82},{1.6,0.67},{2.4,0.55},{3.2,0.45},{4.,0.37},{4.8,0.30},{5.6,0.25},{6.4,0.20},{7.2,0.17},{8.,0.14},{8.8,0.11},{9.6,0.091},{10.4,0.07},{11.2,0.061}};

ListPlot[{golgiLAT,golgiAllL,golgiLATtmd,golgiAllLtmd},Joined->True,PlotTheme->{"Scientific","BoldColor","SansLabels","OpenMarkers"},FrameLabel->{"time [h]","Residence fraction"},PlotLegends->{"LAT","All-L","LAT-TMD","All-L-TMD"}]

dataListT={{erLAT,golgiLAT},{erLATtmd,golgiLATtmd},{erAllL,golgiAllL},{erAllLtmd,golgiAllLtmd}};
dataList=Flatten[dataListT,1];
fitData=Flatten[Table[Table[Map[{j,i,#[[1]],#[[2]]}&,dataListT[[j,i]]],{i,Length[dataListT[[j]]]}],{j,Length[dataListT]}],2];

color={ColorData[106][1],ColorData[106][2]};
modelLabel="ER/Golgi exit";
legendData={"ER data","Golgi data"};
plotLabel={"LAT","LAT-tmd","All-L","All-L-tmd"};
maxTList=Table[SortBy[Flatten[dataListT[[i]],1],First][[-1,1]],{i,Length[dataListT]}];
pdataNoLines=Table[ListPlot[dataListT[[i]],PlotRange->{{-0.02*maxTList[[i]],1.02*maxTList[[i]]},All},PlotStyle->Map[Lighter[#,0.5]&,color],PlotTheme->{"Scientific","BoldColor","SansLabels","OpenMarkers"},FrameLabel->{"time [h]","Residence fraction"},PlotLegends->Placed[legendData,{0.75,0.75}],PlotLabel->plotLabel[[i]],ImageSize->300,AspectRatio->1],{i,Length[dataListT]}];
Magnify[GraphicsGrid[Partition[pdataNoLines,2],ImageSize->{800,800}],1]

Joint analysis

(* list of all model parameters that go into the function, in order *)
modelCur=modelAllJointLogPS4j;
allModelParams=allModelParamsLogPS4j

(* set the fixed parameters and their values *)
fixedParams=Map[{#[[1]],N[#[[2]],numDigitsPrecision]}&,{{logkϵ,Log[10,10^5]},{logkγ,Log[10,10^5]},{logkd,Log[10,10^-5]},{logkm1Full,Log[10,10^5]},{logkm1tmd,Log[10,10^5]},{logk2,Log[10,10^5]},{logkb,Log[10,10^-5]}}];
(* set the adjustable parameters and their initial values *)
adjParamsInit=Map[{#[[1]],N[#[[2]],numDigitsPrecision]}&,{{logKELAT,Log[10,3]},{logKEAllL,Log[10,0.5]},{logKEtmd,Log[10,0.01]},{logKGLAT,Log[10,2]},{logKGAllL,Log[10,0.1]},{logka,Log[10,3]},{logkc,Log[10,2]}}];
setModelParams[]
plotModelPrediction[]

maxIterations=500;
method="LevenbergMarquardt";
c=0;
timeInit=AbsoluteTime[];
nlm=NonlinearModelFit[fitData,Apply[modelCur,Join[nlmParams,{σ,γ,t}]],adjParamsInit,{σ,γ,t},Method->method,MaxIterations->maxIterations,EvaluationMonitor:>{c++}];
timeFinal=AbsoluteTime[];
bfp=nlm["BestFitParameters"];
adjParamNames=Map[#[[1]]&,nlm["BestFitParameters"]];
TableForm[N[Map[{StringTake[ToString[#[[1]]],{4,-1}],10^#[[2]],10^#[[1]]/.bfp}&,adjParamsInit]],TableHeadings->{None,{"param","initial","optimized"}}]
Print["elapsed time: ",ToString[Round[timeFinal-timeInit,0.1]]," s"];
Print["number of function calls: ",ToString[c]];
Print["sum of squared residuals: ",ToString[Total[nlm["FitResiduals"]^2]]];
plotModelFit[]

linearCITable=Quiet[MapThread[{StringTake[ToString[#1[[1]]],{4,-1}],10^#2[[1]],Chop[10^#2[[3,1]]],10^#2[[3,2]]}&,{nlm["BestFitParameters"],nlm["ParameterConfidenceIntervalTableEntries"]}]];
TableForm[linearCITable,TableHeadings->{None,{"Parameter","Best fit (linear)","Min CI","Max CI"}}]

logCITable=Quiet[MapThread[{StringTake[ToString[#1[[1]]],{4,-1}],#2[[1]],Chop[#2[[3,1]]],#2[[3,2]]}&,{nlm["BestFitParameters"],nlm["ParameterConfidenceIntervalTableEntries"],Range[Length[nlm["BestFitParameters"]]]}]];
TableForm[logCITable,TableHeadings->{None,{"Parameter","Best fit (log)","Min CI","Max CI"}}]

corrMat=Round[nlm["CorrelationMatrix"],0.01];
TableForm[corrMat,TableHeadings->ConstantArray[Map[#[[1]]&,nlm["BestFitParameters"]],2]]