% Kippner 2018 - IL-2 model helper code % 2018 May: This version sets the intial value for IL2RA and IL2RB/G to a random value based on flow expression data. % Use with simbio model where these two speciel values are set to 0. % Flow data in file export_Samples_CD25 Cy5 and Tube2 CD1222 FITC_002_014.xlsx % This is a function to control the addition of ligand to the simbio model % pulseTimes = [] defines what time points are used. % L = [] defines the amount of L added/removed at each time point. Here it is 100pM % Ctrl+C stops a run %%% Species in order in model, with default starting value: % 1 Rs 1500 % 2 L 0 % 3 Cs 10 % 4 Ri 300 % 5 Li 0.01 % 6 Ci 1 % 7 Ld 1 % 8 Y 4.22465e+09 % 9 Cs_total 0 % 10 Rs_occ 0 % 11 A 0 (set to 0 for this version) % 12 BG 0 (set to 0 for this version) % 13 LA 1 function f = nextPulseTime100pM2018VarySubunit() clear; % Clear workspace. Note: do not close all; we want to graph all data in same figure Filename = ['ModelResults100pM_RandomizingIL2RAandIL2RBG_30s30s_',datestr(now,30),'.mat']; % Name of results file mdl = sbioloadproject('Kippner2018.sbproj'); %load sbio model (version w IL2RA and IL2RBG values=0 get(mdl.m1.getconfigset('active').SolverOptions) % QC: check what the solver settings are (AbsTol default in 1e-12) mdl.m1.getconfigset('active').SolverOptions.AbsoluteTolerance = 1e-11; % try this instead for 5m5m setting to make sure integration tolerance is met species = get(mdl.m1,'Species'); %define where species are %%% == Generate and assign random starting values to receptor subunits within a range based on experimental data === % Flow data has been normalized 0 to 1 and multiplied by max number [num,text,raw] = xlsread('CD25_Cy5_CD1222_FITC.xlsx','IL2RA and IL2RB'); alpha = num(:,8); % beta = num(:,9); for i = 1:200 rand_id = ceil(rand*length(alpha)); % pick random index between 1 and 10000 (number of data points) alpha_choice = alpha(rand_id); % Find random value for IL2RA within the flow data set beta_choice = beta(rand_id); set (species, 'ConstantAmount', false); set(species(11),'InitialAmount',alpha_choice); % set IL2RA to value defined above set(species(12),'InitialAmount',beta_choice); % set IL2RBG to value defined above %%% From below, the code is correct %%% === Varying the IL2 input. Choose option of interest below === %%% Constant IL2 100pM % pulseTimes = [60]; % L = [1E-10]; %%% 30SEC PULSES pulseTimes = [60 90 120 150 180 210 240 270 300 330 360 390 420 450 480 510 540 570 600 ... 630 660 690 720 750 780 810 840 870 900 930 960 990 1020 1050 1080 1110 1140 1170 1200 ... 1230 1260 1290 1320 1350 1380 1410 1440 1470 1500 1530 1560 1590 1620 1650 1680 1710 1740 1770 1800 ... 1830 1860 1890 1920 1950 1980 2010 2040 2070 2100 2130 2160 2190 2220 2250 2280 2310 2340 2370 2400 ... 2430 2460 2490 2520 2550 2580 2610 2640 2670 2700 2730 2760 2790 2820 2850 2880 2910 2940 2970 3000 ... 3030 3060 3090 3120 3150 3180 3210 3240 3270 3300 3330 3360 3390 3420 3450 3480 3510 3540 3570 3600]; %30s pulse 30s pause L = [1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, ... 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, ... 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, ... 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, ... 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, ... 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10]; %100pM IL2 pulses %%% 1MIN PULSES % pulseTimes = [60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 ... % 1440 1500 1560 1620 1680 1740 1800 1860 1920 1980 2040 2100 2160 2220 2280 2340 2400 2460 2520 2580 2640 ... % 2700 2760 2820 2880 2940 3000 3060 3120 3180 3240 3300 3360 3420 3480 3540 3600]; %time in seconds. 1min on 1min off % L = [10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, ... % 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, ... % 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0, 10e-12, 0]; %%% 2MIN PULSES % pulseTimes = [60 180 210 330 360 480 510 630 660 780 810 930 960 1080 1110 1230 1260 1380 1410 1530 1560 1680 1710 1830 ... % 1860 1980 2010 2130 2160 2280 2310 2430 2460 2580 2610 2730 2760 2880 2910 3030 3060 3180 3210 3330 3360 3480 3510]; %2min pulse 30s pause % L = [1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0 ... % 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0 ... % 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10]; %2m30s %%% 5MIN PULSES % pulseTimes = [60 360 660 960 1260 1560 1860 2160 2460 2760 3060 3360]; %5min pulse 5min pause % L = [1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0, 1E-10, 0]; %5m5m % %%% Run through all the pulses of Il2 specified above for i = 1 : numel(pulseTimes) % loop over all pulse times and add the corresponding L value for each eventTrig=sprintf('time >= %d',pulseTimes(i)); % the event trigger is the time point pulseTimes eventDo=sprintf('L =%d',L(i)); % the event is the addition of L eventObj=addevent(mdl.m1,eventTrig,eventDo); % add event to model at the trigger specified end % %%% Tracking and graphing species of interest simData = sbiosimulate(mdl.m1); % run model with the above info for output simData simData5 = selectbyname(simData,{'Cs'}) % find Cs data in results %%% NOTE: Depending on the model verion used, the species you want may be in a different column than 5. Open the x matrix to check. %%% Having the Cs values from the simbio model, we can now calculate Sr over time % values for logistic delay fcn below y=1; % intercept. Is 1 for normalized Sr b=0.0035; % Cs scaling d=1.8; ku1=0.0025; % steepness ti=2000; % inflection point %%% calculate Sr = STAT5-GFP ratio: Sr=y+(b*simData5.Data)./(d+exp(-ku1*(simData5.Time-ti))); % figure(1) % plot Cs % Don't use this, or there will be a new figure for % every run % hold all; % plot all Cs curves in same figure % sbioplot(simData5); % generate plot of output, choosing from the options defined above % xlabel('Time (s)'); % ylabel('R-L complexes/cell'); figure(2) % plot STAT5 ratio in figure number 3 hold all plot(simData5.Time,Sr) xlabel('Time (s)'); ylabel('Nuclear/Cytosolic STAT5-GFP'); title('STAT5 translocation with pulsatile IL-2 100pM 30s30s for varying IL2RA and IL2BG levels kx0.00175 Vs 0.25 200 runs'); %%% Tracking metrics of interest - Time (min) of max # of Cs per cell and Area under Cs curve x = simData.Data; % create array of output data since the max fcn needs a data array to search in t = simData.Time; % create array of time for above data [Csmax,ind] = max(x(:,3)); % find max Cs (Cs is species in column 3, see simData>DataNames or the list in the model) Timemax = t(ind); % find the time when Cs is at max Int = trapz(t,x(:,1)); % calculate area under Cs curve. Trapz fcn estimates area as trapeziods. [STAT5max, ind2] = max(Sr); % find max for STAT5 ratio Timemax2 = t(ind2); % find time at max STAT5 ratio Int2 = trapz(t,Sr); % calculate area under STAT5 ratio curve % Optional QC: print values of interest for run % Csmax % Int % Timemax % STAT5max % Int2 % Timemax2 end %save(Filename) % Save results with the name specified at the beginning of this file end %%% Troubleshooting help: %equations = getequations(mdl.m1); %Get the equations that Matlab generates for the model.