Added all files

.gitignore

+*.dat
+*.asv
+*.autosave
+*.xlsx
+*.xls
+*.png
+*.fig
+*.pdf
+*.h
+*.c
+*.slxc
+*.out
+*.log
+Parameters/ESS*
+Parameters/PL*
+Plots/*
+Data/SampledErrors_*/*
+Resultfigures*
\ No newline at end of file

Data/addErrosToRealData.m

+clear
+close all
+
+% Load and format data for 80 subjects, and add the estimated combined errors.
+for p = 1:80
+    % Load the clinical data
+    load(['inputdata_avatarmodel_P' num2str(p) '.mat'])
+
+    % load true data, which contains the combined errors
+    load('dataSimulated.mat','dataSimulated')
+    trueData = dataSimulated.simRandomSystematic;
+
+    try
+        %add mean values & time to the struct containing the final
+        %estimation data
+        estimationData = struct;
+        estimationData.MV.mean = data.MV(:,2);
+        estimationData.MV.time = data.MV(:,1);
+        estimationData.AV.mean = data.AV(:,2);
+        estimationData.AV.time = data.AV(:,1);
+        estimationData.AA.mean = data.AC(:,2);
+        estimationData.AA.time = data.AC(:,1);
+        estimationData.PV.mean = data.PV(:,2);
+        estimationData.PV.time = data.PV(:,1);
+
+        estimationData.indtdiast = data.indtdiast;
+        estimationData.extra.tdiast = paramdata.tdiast;
+        estimationData.extra.indtdiast = paramdata.indtdiast;
+
+
+        %add estimated combined data uncertainty
+        [e.estimatedRRerror] = calcRRsmoothingerror(estimationData); %use empirical function to calculate the RR error
+        estimationData.MV.sem = trueData.MV.eRand + trueData.MV.eSyst + abs(e.estimatedRRerror.MV.mean);
+        estimationData.AV.sem = trueData.AV.eRand + trueData.AV.eSyst + abs(e.estimatedRRerror.AV.mean);
+        estimationData.AA.sem = trueData.AA.eRand + trueData.AA.eSyst + abs(e.estimatedRRerror.AA.mean);
+        estimationData.PV.sem = trueData.PV.eRand + trueData.PV.eSyst + abs(e.estimatedRRerror.PV.mean);
+
+
+        estimationData.SBP.sem = trueData.SBP.sem;%inherit the sem
+        estimationData.DBP.sem = trueData.DBP.sem;%inherit the sem
+        estimationData.extra.ESV.sem = trueData.extra.ESV.eRand;
+
+        estimationData.SBP.mean = paramdata.SBP;
+        estimationData.DBP.mean = paramdata.DBP;
+        estimationData.extra.ESV.mean = paramdata.ESV_seg;
+
+        estimationData.extra.EDV = NaN;
+        estimationData.extra.LaESV = NaN;
+
+        % information for simulating the data
+        estimationData.meta.simulationFunction = trueData.meta.simulationFunction;
+        estimationData.meta.modelName = trueData.meta.modelName;
+        estimationData.meta.paramsToOptimize = trueData.meta.paramsToOptimize;
+
+        % Add sem for all data-based parameters
+        estimationData.extra.T.mean = paramdata.T;
+        estimationData.parameters.Caa.mean = paramdata.Caa;
+        estimationData.parameters.ELCo.mean = paramdata.ElCo;
+
+        [estimationData.parameters.Ctot.mean,estimationData.parameters.Rtot.mean,estimationData.parameters.Emax_LV.mean] = calculateDatabasedParameters(estimationData);
+
+        pnames = fieldnames(trueData.parameters);
+        for i = 1:length(pnames)
+            estimationData.parameters.(pnames{i}).sem = trueData.parameters.(pnames{i}).sem;%inherit the sem
+        end
+
+        %% Optional: create plot of the resulting estimation data to make sure it is reasonable
+        figure()
+        nexttile
+        title('mv')
+        errorbar(estimationData.MV.time,estimationData.MV.mean,estimationData.MV.sem)
+        nexttile
+        title('av')
+        errorbar(estimationData.AV.time,estimationData.AV.mean,estimationData.AV.sem)
+        nexttile
+        title('aa')
+        errorbar(estimationData.AA.time,estimationData.AA.mean,estimationData.AA.sem)
+        nexttile
+        title('pv')
+        errorbar(estimationData.PV.time,estimationData.PV.mean,estimationData.PV.sem)
+        nexttile
+        title('parameters')
+        hold on
+        for i = 1:length(pnames)
+            errorbar(i,estimationData.parameters.(pnames{i}).mean,estimationData.parameters.(pnames{i}).sem,'*');
+        end
+
+        nexttile
+        hold on
+        errorbar(1,estimationData.SBP.mean,estimationData.SBP.sem)
+        errorbar(2,estimationData.DBP.mean,estimationData.DBP.sem)
+        title('bp')
+
+        %% Save the created estimation data for this subject
+        save(['dataP' num2str(p) '.mat'],'estimationData')
+
+    catch
+        disp(['couldnt do p ' num2str(p)])
+    end
+
+end
+
+%% Plot selected data
+load('dataP78.mat','estimationData')
+cdata = estimationData;
+load('dataP33.mat','estimationData')
+pdata = estimationData;
+
+figure()
+nexttile
+title('mv')
+hold on
+errorbar(cdata.MV.time./cdata.MV.time(end),cdata.MV.mean,cdata.MV.sem)
+errorbar(pdata.MV.time./pdata.MV.time(end),pdata.MV.mean,pdata.MV.sem)
+
+nexttile
+title('av')
+hold on
+errorbar(cdata.AV.time./cdata.MV.time(end),cdata.AV.mean,cdata.AV.sem)
+errorbar(pdata.AV.time./pdata.MV.time(end),pdata.AV.mean,pdata.AV.sem)
+
+nexttile
+hold on
+title('aa')
+errorbar(cdata.AA.time./cdata.MV.time(end),cdata.AA.mean,cdata.AA.sem)
+errorbar(pdata.AA.time./pdata.MV.time(end),pdata.AA.mean,pdata.AA.sem)
+
+nexttile
+hold on
+title('pv')
+errorbar(cdata.PV.time./cdata.MV.time(end),cdata.PV.mean,cdata.PV.sem)
+errorbar(pdata.PV.time./pdata.MV.time(end),pdata.PV.mean,pdata.PV.sem)
+
+nexttile
+title('parameters')
+hold on
+for i = 1:length(pnames)
+    errorbar(i,cdata.parameters.(pnames{i}).mean,cdata.parameters.(pnames{i}).sem,'*');
+    errorbar(i+0.2,pdata.parameters.(pnames{i}).mean,pdata.parameters.(pnames{i}).sem,'*');
+end
+xticks(1:length(pnames))
+xticklabels(pnames)
+
+nexttile
+hold on
+errorbar(1,cdata.SBP.mean,cdata.SBP.sem)
+errorbar(2,cdata.DBP.mean,cdata.DBP.sem)
+errorbar(1.2,pdata.SBP.mean,pdata.SBP.sem)
+errorbar(2.2,pdata.DBP.mean,pdata.DBP.sem)
+title('bp')

Data/calcRRsmoothingerror.m

+function [error] = calcRRsmoothingerror(data)
+% Empirical function to estimate the smoothing error due to RR variability
+% in 4D flow MRI data. The error is calculated based on the blood flow curves 
+% in MV, PV, AV and AA in the given dataset. This is used to estimate data
+% uncertainty when not knowing the true error.
+% This function is based on 40 data points per flow curve.
+
+datacontent = fieldnames(data);
+error = struct;
+%% MV
+if ismember('MV',datacontent)
+    tstep = max(1,round(length(data.MV.time)/40));
+    diastdata = data.MV.mean(data.indtdiast:end);
+    [~,peakinds] = findpeaks(diastdata,'NPeaks',2,'MinPeakProminence',2,'MinPeakDistance',4);
+    peakinds = peakinds+data.indtdiast-1;
+    calculatedDiffMV = zeros(size(data.MV.mean));
+    calculatedDiffMV(data.indtdiast:end) = -0.10.*data.MV.mean(data.indtdiast:end);
+
+    % The E and A waves are affected the most
+    if ~isempty(peakinds) %E-wave is affected more than the A-wave (negatively)
+        calculatedDiffMV(peakinds(1)-2*tstep:peakinds(1)+3*tstep) = -0.15.*data.MV.mean(peakinds(1)-2*tstep:peakinds(1)+3*tstep);
+    end
+    if length(peakinds) > 1 % A-wave (positively)
+        if peakinds(2) < 38
+            calculatedDiffMV(peakinds(2)-2*tstep:peakinds(2)+3*tstep) = 0.13.*data.MV.mean(peakinds(2)-2*tstep:peakinds(2)+3*tstep);
+        else
+            disp('A peak too close to timevector end, no extra sd added')
+        end
+    end
+
+    % The start of diastole (pos)
+    calculatedDiffMV(data.indtdiast:data.indtdiast+5*tstep) = 15;
+
+    error.MV.mean = calculatedDiffMV;
+    error.MV.time =  data.MV.time;
+end
+
+%% PV
+%only a small smoothing error, approximated with +-5 ml/s
+if ismember('PV',datacontent)
+    error.PV.mean = 5.*ones(size(data.PV.mean));
+    error.PV.time = data.PV.time;
+
+    %find the PV peaks (for a healthy "textbook" flow)
+    [pks,locs] = findpeaks(data.PV.mean,'Npeaks',2);
+    Sind = locs(1);
+    Dind = locs(2);
+    [pks,locs] = findpeaks(-data.PV.mean,'Npeaks',2);
+    Aind = locs(2);
+    midSD = locs(1);
+
+    tstep = max(1,round(length(data.MV.time)/40));
+
+    %positive error around start and in the end after the A wave (leave as is)
+    
+    % negative error around S, D, A
+    error.PV.mean(Sind-2*tstep:Sind+4*tstep) =  -error.PV.mean(Sind-2*tstep:Sind+4*tstep);
+    error.PV.mean(Dind-2*tstep:Aind+2*tstep) =  -error.PV.mean(Dind-2*tstep:Aind+2*tstep);
+
+    % positive error in between S % D (leave as is)
+end
+
+%% AV
+if ismember('AV',datacontent)
+    tstep = max(1,round(length(data.AV.time)/40));
+    calculatedDiffAV = zeros(size(data.AV.mean));
+    calculatedDiffAV(1:data.indtdiast) = -0.03.*data.AV.mean(1:data.indtdiast);% -3% in general, more errors at peak
+    calculatedDiffAV(data.indtdiast:data.indtdiast+2*tstep) = 10;%positive error
+    error.AV.mean = calculatedDiffAV;
+    error.AV.time =  data.AV.time;
+end
+
+%% AA
+if ismember('AA',datacontent)
+    tstep = max(1,round(length(data.AA.time)/40));
+    calculatedDiffAA = zeros(size(data.AA.mean));
+    calculatedDiffAA(1:data.indtdiast) = -0.03.*data.AA.mean(1:data.indtdiast);%-3 in general
+    calculatedDiffAA(data.indtdiast:data.indtdiast+2*tstep) = 10;
+
+    error.AA.mean = calculatedDiffAA;
+    error.AA.time =  data.AA.time;
+end
+
+end
\ No newline at end of file

Data/calculateDatabasedParameters.m

+function [Ctot,Rtot,Emax_LV] = calculateDatabasedParameters(data)
+% Calculate parameters on a subject-specific basis based on the input data
+
+%% Emax LV
+rho_blood=1.06; %g/ml (constant)
+B_av=rho_blood/(2*data.parameters.ELCo.mean^2);
+maxQav = max(data.AV.mean);
+deltaPconst = 0.06/133.322;
+deltaP_av = B_av*(maxQav^2)*deltaPconst; %Pressure gradient accross the aortic valve
+k=0.9; %Constant to relate systolic blood pressure and end systolic blood pressure: Pes=k*Ps 
+V0_LV=10; %(mL) Mynard et al. (constant)
+Emax_LV=k*(data.SBP.mean+deltaP_av)/(data.extra.ESV.mean-V0_LV); %Maximal elastance of the LV
+
+
+%% Ctot
+PP=data.SBP.mean-data.DBP.mean;    
+meanQav = mean(data.AV.mean);
+Ctot=meanQav/PP;  %Total compliance of the system
+
+
+%% Rtot
+HR=60*1/data.extra.T.mean;  %Heart rate (beats/min)
+MAP=data.DBP.mean+(1/3+HR*0.0012)*(data.SBP.mean-data.DBP.mean); % Beléns calculation of MAP, takes HR in consideration azminia M, Trivedi A, Molnar J, Elbzour M, Guerrero M, Salem Y, Ahmed A, Khosla S, Lubell DL . Validation of a new formula for mean arterial pressure calculation: the new formula is superior to the standard formula. Catheter Cardiovasc Interv 2004; 63: 419–425.
+Rtot=MAP/meanQav;     %Total resistance of the system
+
+
+end
\ No newline at end of file

Data/createSampledMeasurementError.m

+function createSampledMeasurementError(numberOfSamples,experimentName)
+% experimentName = 'simRandomSystematic' or  simRandomOnly
+rng("default") %for reproducibility
+
+basefolder = split(pwd,'Uncertainty-estimation');
+basefolder = fullfile(basefolder{1},'Uncertainty-estimation');
+addpath(genpath(basefolder))
+
+savefolder = fullfile(basefolder,'Data',['SampledErrors_' experimentName]);
+mkdir(savefolder)
+
+load('dataSimulated.mat','dataSimulated')
+trueData = dataSimulated.(experimentName);
+
+for n = 1:numberOfSamples
+    disp(100*n/numberOfSamples)
+    [~,e] = sampleData(trueData);
+    % [~,e] = sampleData(trueData,[],n); %load RR error
+    save(fullfile(savefolder,sprintf('sampledMeasurementError%d.mat',n)),'e')
+end
+
+fprintf('%d samples created and saved in %s\n',numberOfSamples,savefolder)
+end
\ No newline at end of file

Data/createTrueData_simulated.m

+%% createTrueData_simulated
+clear
+% Script to create estimation data based on a model simulation as "true" data
+% also add the distribution of errors, that can be sampled using
+% createSampledMeasurementError (in sampleData)
+
+% NOTE: data names cannot contain underscore!
+
+% add paths
+basefolder = split(pwd,'Uncertainty-estimation');
+basefolder = fullfile(basefolder{1},'Uncertainty-estimation');
+
+%add all project folders to the matlab path
+addpath(genpath(basefolder))
+
+% setup AMICI toolbox
+run(fullfile(basefolder, 'Requirements', 'AMICI-0.10.11_SS_eventFix', 'matlab', 'installAMICI.m'))
+
+
+%% Create the "true" simulation data
+
+%setup
+% Simulation settings
+% sensi: sensitivity order - 0 = no calculation of sensitivities (default)
+% maxsteps: maximum number of integration steps
+simulationoptions = amioption('sensi',0,'maxsteps',1e4);
+% set simulation tolerances
+simulationoptions.atol = 1e-16;
+simulationoptions.rtol = 1e-8;
+
+load('trueParameters.mat','trueParams','trueConstants','truedata','ind','paramdata','extradata')
+paramNamesOld = {'Cpvc' 'Rpu' 'Rpv' 'Lpv' 'Rtot' 'Ctot' 'ElCo' 'Caa' 'Emax_LA' 'Emax_LV' 'Emin_LA' 'Emin_LV' 'Lao' 'Lav'...
+    'Lmv' 'Ppu' 'Rao' 'Rmv' 'k_diast_LA' 'k_diast_LV' 'k_syst_LA'...
+    'k_syst_LV' 'm1_LA' 'm1_LV' 'm2_LA' 'm2_LV' 'onset_LA' 'onset_LV'};
+constantsNamesOld = {'aaCorr' 'avCorr' 'mvCorr' 'tdiast' 'Ks_LA' 'Ks_LV'...
+    'V0_LA' 'V0_LV' 'RLAvisc' 'RLVvisc' 'Raa' 'Rpc' 'Rpvc' 'T' 'rho_blood' 'norm_factor_LA' 'norm_factor_LV'};
+
+%load parameter names and indexes
+T = trueConstants(ind.T);
+simulatedDataTrue.extra.T.mean = T;
+simulatedDataTrue.extra.indtdiast = paramdata.indtdiast;
+simulatedDataTrue.extra.tdiast = paramdata.tdiast;
+simulatedDataTrue.parameters.ELCo.mean = trueParams(ind.ElCo);
+simulatedDataTrue.parameters.Caa.mean = trueParams(ind.Caa);
+simulatedDataTrue.parameters.Ctot.mean = trueParams(ind.Ctot);
+simulatedDataTrue.parameters.Rtot.mean = trueParams(ind.Rtot);
+simulatedDataTrue.parameters.Emax_LV.mean = trueParams(ind.Emax_LV);
+oldind = ind;
+[~,~,paramNames,constantsNames,units,ind] = loadParameters(simulatedDataTrue);
+
+%convert the true params to the new order or params and constants
+trueConstants = trueConstants(4:end); %remove corr constants
+trueParams = [trueParams,40,40,40,40]; % add correction parameters
+
+%simulate
+simulationoptions.x0 = truedata.IC;
+datatime = truedata.time;
+simtime = sort([datatime,0:0.001:T]);
+simtime = unique(simtime);
+datatimeinds = find(ismember(simtime,datatime));
+cd ../Simulation
+warning('off')
+[sol] = simulate_avatar_ss(trueParams,trueConstants,simulationoptions,ind,simtime,'avatar_corr');
+warning('on')
+cd ../Data
+%set the simulated data based on the simulation
+simulatedDataTrue.MV.time = sol.t;
+simulatedDataTrue.AV.time = sol.t;
+simulatedDataTrue.AA.time = sol.t;
+simulatedDataTrue.PV.time = sol.t;
+simulatedDataTrue.MV.mean = sol.x(:,ind.MV);
+simulatedDataTrue.AV.mean = sol.x(:,ind.AV);
+simulatedDataTrue.AA.mean = sol.x(:,ind.AA);
+simulatedDataTrue.PV.mean = sol.x(:,ind.PV);
+simulatedDataTrue.extra.SV_MV.mean = trapz(sol.t,sol.x(:,ind.MV));
+simulatedDataTrue.extra.SV_AV.mean = trapz(sol.t,sol.x(:,ind.AV));
+simulatedDataTrue.extra.SV_AA.mean = trapz(sol.t,sol.x(:,ind.AA));
+simulatedDataTrue.extra.EDV.mean = max(sol.x(:,ind.LV));
+simulatedDataTrue.extra.ESV.mean = min(sol.x(:,ind.LV));
+simulatedDataTrue.extra.LaESV = max(sol.y(:,ind.Vla));
+simulatedDataTrue.Vlv = sol.x(:,ind.LV);
+d.SBP.mean = truedata.SBP;
+[simulatedDataTrue.SBP.mean,simulatedDataTrue.DBP.mean] = brachialpressure(sol,ind,d);
+simulatedDataTrue.SBP.time = NaN;simulatedDataTrue.DBP.time = NaN;
+
+%% Calculate RR smoothing
+rrNumber = 400;
+randomDistribution = 0;
+fprintf('Simulating %d heartbeats...\n',rrNumber)
+[sol,solTrue,~] = simulate_avatar_RR(trueParams,trueConstants,simulationoptions,ind,simtime,'avatar_corr',rrNumber,randomDistribution,0); 
+rrdata.MV.mean = sol.x(:,ind.MV);
+rrdata.AV.mean = sol.x(:,ind.AV);
+rrdata.AA.mean = sol.x(:,ind.AA);
+rrdata.PV.mean = sol.x(:,ind.PV);
+rrdata.MV.time = sol.t;
+rrdata.AV.time = sol.t;
+rrdata.AA.time = sol.t;
+rrdata.PV.time = sol.t;
+[~,rrdata.indtdiast] = min(abs(sol.t-paramdata.tdiast));
+[estimatedRRerror] = calcRRsmoothingerror(rrdata); %re-estimated for each sampled data in SampleData (createSampledMeasurementError)
+
+%% Add uncertainty: only random error included in sigma (sem)
+exp = simulatedDataTrue;
+exp.meta.timeResolution = 40;
+exp = setTemporalResolution(exp,exp.meta.timeResolution);
+
+% Blood flow curves: standard deviation of the random error
+exp.MV.eRand = abs(exp.MV.mean).*0.0539; % 5.4% of flow  - different in each timepoint
+exp.MV.eSyst= 0.9*7; % 0.9 cm/s * vessel area  (doi:10.1016/j.jsha.2010.07.007, https://doi.org/10.1152/ajpheart.00269.2004 ) - same for all timepoints
+exp.AV.eRand = abs(exp.AV.mean).*0.0539; % 5.4% of flow
+exp.AV.eSyst= 0.9*5; % 0.9 cm/s * vessel area (doi: 10.4137/CMC.S15716, doi:10.1053/euhj.2001.2782,  pi*(3/2)^2)
+exp.AA.eRand = abs(exp.AA.mean).*0.0539; % 5.4% of flow
+exp.AA.eSyst= 0.9*7;% 0.9 cm/s * vessel area (doi:10.1053/euhj.2001.2782   pi*(2.5/2)^2)
+exp.PV.eRand = abs(exp.PV.mean).*0.0806; % 12% of flow
+exp.PV.eSyst = 0.9*5.3;% 0.9 cm/s * vessel area (4*pi*(1.3/2)^2 = 5.3 cm^2  (Wittkampf et al., 2003; Kim et al., 2005))
+
+% where the flow is close to 0, the % error is too small. Correct this:
+exp.MV.eRand(exp.MV.eRand < mean(exp.MV.eRand)) = mean(exp.MV.eRand);
+exp.AV.eRand(exp.AV.eRand < mean(exp.AV.eRand)) = mean(exp.AV.eRand);
+exp.AA.eRand(exp.AA.eRand < mean(exp.AA.eRand)) = mean(exp.AA.eRand);
+exp.PV.eRand(exp.PV.eRand < mean(exp.PV.eRand)) = mean(exp.PV.eRand); 
+
+% add the error to sigma (only the random error)
+exp.MV.sem = exp.MV.eRand; 
+exp.AV.sem = exp.AV.eRand;
+exp.AA.sem = exp.AA.eRand;
+exp.PV.sem = exp.PV.eRand;
+
+% Data-based parameters
+exp.parameters.Caa.eRand = simulatedDataTrue.parameters.Caa.mean.*0.15; % 15% error
+exp.parameters.ELCo.eRand = simulatedDataTrue.parameters.ELCo.mean.*0.12; % 12% error
+exp.parameters.Ctot.eRand = simulatedDataTrue.parameters.Ctot.mean.*0.20; % 25%, calculated from other variables
+exp.parameters.Rtot.eRand = simulatedDataTrue.parameters.Rtot.mean.*0.082; % 12%, calculated from other variables
+exp.parameters.Emax_LV.eRand = simulatedDataTrue.parameters.Emax_LV.mean.*0.077; % 10%, calculated from other variables
+
+exp.parameters.Caa.eSyst = 0; % no systematic errors
+exp.parameters.ELCo.eSyst = 0; % no systematic errors
+exp.parameters.Ctot.eSyst = 0; % no systematic errors added, the parameter is calculated from other measurements
+exp.parameters.Rtot.eSyst = 0; % no systematic errors added, the parameter is calculated from other measurements
+exp.parameters.Emax_LV.eSyst = 0; % no systematic errors added, the parameter is calculated from other measurements
+
+exp.parameters.Caa.sem = exp.parameters.Caa.eRand;
+exp.parameters.ELCo.sem = exp.parameters.ELCo.eRand;
+exp.parameters.Ctot.sem = exp.parameters.Ctot.eRand;
+exp.parameters.Rtot.sem = exp.parameters.Rtot.eRand;
+exp.parameters.Emax_LV.sem = exp.parameters.Emax_LV.eRand;
+
+% Brachial pressure
+exp.SBP.eRand = 8; % 8 mmHg error
+exp.SBP.eSyst= 0;
+exp.SBP.sem = exp.SBP.eRand;
+
+exp.DBP.eRand = 8; % 8 mmHg error
+exp.DBP.eSyst = 0;
+exp.DBP.sem = exp.DBP.eRand;
+
+% LV segmentation
+exp.extra.ESV.eRand = simulatedDataTrue.extra.ESV.mean * 0.05; % 5% random error
+exp.extra.ESV.eSyst= 0;
+exp.extra.ESV.sem = exp.extra.ESV.eRand;
+
+% Other information
+exp.meta.RRnum = 400; %average of 400 random heartbeats
+exp.meta.addRRtoSem = 0; %do not add the systematic RR error to sigme
+exp.meta.description = 'only random errors in sigma';
+exp.meta.paramsToOptimize = 1:length(paramNames);%include the corr parameters
+exp.meta.simulationFunction = 'simulate_avatar_correction';
+exp.meta.modelName = 'avatar_corr';
+expname = 'simRandomOnly';
+dataSimulated.(expname) = exp;
+
+%% Add uncertainty: both random and systematic errors in the sigma/sem
+%calculate rr error
+sind = find(ismember(solTrue.t,exp.MV.time));
+[~,sindunique] = unique(solTrue.t(sind));
+sind = sind(sindunique);
+datanames = {'MV','AV','AA','PV'};
+for n = 1:length(datanames)
+     estimatedRRerror.(datanames{n}).mean = estimatedRRerror.(datanames{n}).mean(sind)';
+end
+
+%all errors are the same as above, it is just the sigma and the simulation
+%function that changes (plus which parameters that are estimated)
+exp.MV.sem = exp.MV.eRand + exp.MV.eSyst + abs(estimatedRRerror.MV.mean');
+exp.AV.sem = exp.AV.eRand + exp.AV.eSyst + abs(estimatedRRerror.AV.mean');
+exp.AA.sem = exp.AA.eRand + exp.AA.eSyst + abs(estimatedRRerror.AA.mean');
+exp.PV.sem = exp.PV.eRand + exp.PV.eSyst + abs(estimatedRRerror.PV.mean');
+
+exp.meta.simulationFunction = 'simulate_avatar_ss';
+exp.meta.modelName = 'avatar_corr';
+exp.meta.paramsToOptimize = 1:find(strcmp(paramNames,'onset_LV'));%do not include the corr parameters
+exp.meta.addRRtoSem = 1; %DO add the systematic RR error to sigma
+expname = 'simRandomSystematic';
+dataSimulated.(expname) = exp;
+
+
+%% Save the full data structure
+% add extra information
+simulatedDataTrue.allSimulations = sol;
+simulatedDataTrue.allParameters = trueParams;
+simulatedDataTrue.constants = trueConstants;
+simulatedDataTrue.ind = ind;
+simulatedDataTrue.options = simulationoptions;
+simulatedDataTrue.simtime = simtime;
+simulatedDataTrue.simFunc = 'simulate_avatar_corr';
+simulatedDataTrue.modelName = 'avatar_corr';
+
+save('dataSimulated.mat','dataSimulated','simulatedDataTrue')
+