3 % model of the human follicle growth
5 % two chemical components are assumed to be characteristic
for the
6 % growth of human follicles.
8 % A deterministic model
for the two species is given by a positive
9 % feedback loop, a bistable
switch, i.e. 2-d nonlinear ODE.
11 % An approximation to the Chemical Master Equation gives a linear
14 % here, the state vector is the probability of states (pair of number of two
15 % chemical components) and the matrix A describes the change in
16 % these probabilities over time. The number of states is
17 % p = (range+1)*(range+1) hence rectangular domain
19 % The output is given by params.output_type:
20 %
'separatrix': an average over a subset of states, the ones below
21 % the diagonal
"separatrix", which are assumed to be still
"in rest",
22 % i.e. waiting
for their growth specified by model.output_range
23 %
'expectation_M': the expectation of M state is computed
24 %
'expectation_N': the expectation of N state is computed
26 % Q: D.h. die Zellen oberhalb der separatrix sind einfach abgestorben, oder zur
29 % Explanations from S. Waldherr
"A stochastic switch for the
30 % primordial to primary transition in follicle growth", IST
31 %
internal report (2009-04-09) v1
33 % Bernard Haasdonk 22.9.2009
38 model = lin_ds_model_default;
39 model.data_const_in_time = 1;
43 %%%%%%%%%%%%%%%%%%%%%%% parameters from steffen
48 if isfield(params,'range')
49 model.range = params.range;
51 model.range = 150; % connected with end time, increase e.g. to 50 later
53 if isfield(params,'output_range')
54 model.output_range = params.output_range;
56 model.output_range = 50; % roughly between two bulbs
58 if (model.output_range > model.range)
59 error('please choose output range <= range');
61 % determines the initial distribution in a lower left triangle
62 if isfield(params,'init_range')
63 model.init_range = params.init_range;
65 model.init_range = 20;
68 if isfield(params,'reg_epsilon')
69 model.reg_epsilon = params.reg_epsilon;
71 model.reg_epsilon = 0;
73 if ~isfield(params,'output_type')
74 params.output_type = 'separatrix';
77 % u1 = 0.01+(i-6)*0.001;
78 % u2 = 0.01+(i-6)*0.0005;
79 % k1 = 4*u1+(i-6)*0.0001;
82 % parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 %model_data = my_ds_gen_model_data(model);
94 % used as components later:
95 %model.k1 = model.k1fact * model.u1;
96 %model.V1 = model.V1fact * model.u1;
97 %model.V2 = model.V2fact * model.u2;
99 model.mu_names = {
'u1',
'u2'};
100 %model.mu_names = {
'u1',
'u2',
'k1fact',
'V1fact',
'V2fact'};
101 model.mu_ranges = {[0.005,0.02],[0.005,0.02]};
103 %model.T = 1000000; % Minutes, i.e. 1.9 Jahre
105 %model.T = 1000000; % Minutes, i.e. 1.9 Jahre
107 model.T = 10000000; % Minutes, i.e. 19 Jahre
109 %model.T = 5000000; % Minutes, i.e. 9.5 Jahre
111 %model.T = 50000000; % Minutes, i.e. 95 Jahre
113 %disp(
'nt to be adjusted later!!!') => 500 OK!!
114 model.affinely_decomposed = 1;
115 model.A_function_ptr = @A_func;
116 % ratio of first to second ellipses axes
117 %model.A_axes_ratio = 2;
118 %model.A_axes_ratio = 2;
119 model.B_function_ptr = @B_func;
120 switch params.output_type
122 model.C_function_ptr = @C_func_separatrix;
124 model.C_function_ptr = @C_func_expectation_M;
126 model.C_function_ptr = @C_func_expectation_N;
128 error(
'params.output_type not known');
130 model.D_function_ptr = @D_func;
131 model.x0_function_ptr = @x0_func;
132 % shift of x0 in third coordinate
134 model.u_function_ptr = @(model) 0; % define u constant
135 %model.u_function_ptr = @(model) sin(4*model.t);
136 %model.u_function_ptr = @(model) 0; % define u constant 0
137 %model.u_function_ptr = @(model) model.t; % define u = 1 * t
138 model.G_matrix_function_ptr = ...
139 @(model,model_data) speye(model.dim_x,model.dim_x);
140 % the following lead to both larger C1 and C2!!!
141 %model.G_matrix_function_ptr = @(model) diag([2,0.5,1]);
142 %model.G_matrix_function_ptr = @(model) diag([0.5,2,1]);
144 % set
function pointers to main model methods
145 model.gen_model_data = @my_ds_gen_model_data;
146 %model.detailed_simulation = @lin_ds_detailed_simulation;
147 model.gen_detailed_data = @my_ds_gen_detailed_data;
148 %model.gen_reduced_data = @lin_ds_gen_reduced_data;
149 %model.rb_simulation = @lin_ds_rb_simulation;
150 %model.plot_sim_data = @my_ds_plot_sim_data;
151 model.plot_sim_data_state = @plot_state_probabilities;
152 model.plot_slice = @my_plot_slice;
153 model.axis_tight = 1;
155 %model.rb_reconstruction = @lin_ds_rb_reconstruction;
156 % determin model data, i.e. matrix components, matrix G
158 model.orthonormalize = @model_orthonormalize_qr;
160 %model_data = gen_model_data(model); % matrix components
161 %model.dim_x = size(model_data.X,2);
162 model.dim_x = (model.range+1)*(model.range+1);
164 %model.decomp_mode = 1;
165 %tmp = model.C_function_ptr(model,[]);
167 %model.dim_y = size(tmp{1},1);
168 %model.decomp_mode = 0;
170 %model.dim_y = size(model_data.C_components{1},1);
172 if isfield(params,
'output_plot_indices')
173 model.output_plot_indices = params.output_plot_indices;
175 model.output_plot_indices = 1:model.dim_y;
178 model.theta = 1; % implicit discretization
180 % be sure to determine correct constants once, if anything in the
181 % problem setting has changed!!
183 %model.error_estimation = 1; % model is capable of error_estimation
184 model.enable_error_estimator = 1;
185 if ~isfield(params,'estimate_bounds')
186 params.estimate_bounds = 0;
188 % the following only needs to be performed, if the model is changed.
189 if (params.estimate_bounds == 0)
191 switch params.output_type
193 model.state_bound_constant_C1 = 1 ;
194 model.output_bound_constant_C2 = sqrt(model.dim_x);
196 model.state_bound_constant_C1 = 1 ;
198 model.output_bound_constant_C2 = 100000;
200 warning('please adjust the computation of the error estimator constants!')
201 model.state_bound_constant_C1 = 100000 ;
202 model.output_bound_constant_C2 = 100000;
204 %elseif params.estimate_bounds == 1
205 % model.estimate_lin_ds_nmu = 3;
206 % model.estimate_lin_ds_nX = 10;
207 % model.estimate_lin_ds_nt = 10;
209 % [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
210 % model.state_bound_constant_C1 = C1;
211 % model.output_bound_constant_C2 = C2;
213 model.estimate_lin_ds_nmu = 3;
214 model.estimate_lin_ds_nX = 10;
215 model.estimate_lin_ds_nt = 10;
217 [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
218 model.state_bound_constant_C1 = C1;
219 model.output_bound_constant_C2 = C2;
222 model.inner_product_matrix_algorithm =@(model,model_data) ...
223 speye(model.dim_x,model.dim_x);
225 model.error_algorithm = @(u1,u2,dummy1,dummy2) ...
226 sqrt(max(sum((u1-u2).^2),0));
229 %model.error_norm = 'l2';
230 %% 'RB_extension_max_error_snapshot'};
231 model.RB_stop_timeout = 12*60*60; % 4 hours
232 model.RB_stop_epsilon = 1e-16;
233 model.RB_error_indicator = 'projection-error'; % true projection error
234 %model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
236 model.RB_stop_Nmax = 200;
237 model.RB_generation_mode = 'greedy_uniform_fixed';
238 model.RB_numintervals = [20,20];
239 model.RB_detailed_train_savepath = ...
240 'follicle_rect_model_detailed_train';
241 model.RB_save_temp_detailed_data = 1;
242 model.orthonormalize = @model_orthonormalize_qr;
245 model.transformed_colormap = @my_transformed_colormap;
247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
248 % functions defining the dynamical system and input/initial data
249 function A = A_func(model,model_data);
250 A = eval_affine_decomp(@A_components,@A_coefficients,...
253 function Acomp = A_components(model,model_data)
259 A1 = sparse([],[],[],p,p,7*p); % for k1
260 A2 = sparse([],[],[],p,p,7*p); % for u1
261 A3 = sparse([],[],[],p,p,7*p); % for u2
262 A4 = sparse([],[],[],p,p,7*p); % for V1
263 A5 = sparse([],[],[],p,p,7*p); % for V2
264 A6 = -speye(p,p); % for regularization
267 j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
270 j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
273 j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
275 A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
276 j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
277 A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
278 A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
279 j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
280 A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
283 Acomp = {A1, A2, A3, A4, A5, A6};
285 function Acoeff = A_coefficients(model)
286 %Acoeff = [model.k1, model.u1, model.u2, model.V1, model.V2];
287 k1 = model.k1fact * model.u1;
288 V1 = model.V1fact * model.u1;
289 V2 = model.V2fact * model.u2;
290 Acoeff = [k1, model.u1, model.u2, V1, V2, model.reg_epsilon];
292 function B = B_func(model,model_data)
293 B = eval_affine_decomp(@B_components,@B_coefficients,model,model_data);
295 function Bcomp = B_components(model,model_data)
296 Bcomp = {zeros(model.dim_x,1)};
298 function Bcoeff = B_coefficients(model)
301 function C = C_func_separatrix(model,model_data)
302 C = eval_affine_decomp(@C_components_separatrix,...
303 @C_coefficients,model,model_data);
305 function C = C_func_expectation_M(model,model_data)
306 C = eval_affine_decomp(@C_components_expectation_M,...
307 @C_coefficients,model,model_data);
309 function C = C_func_expectation_N(model,model_data)
310 C = eval_affine_decomp(@C_components_expectation_N,...
311 @C_coefficients,model,model_data);
313 function Ccomp = C_components_separatrix(model,model_data)
314 i = find((model_data.Ms+model_data.Ns)<=model.output_range);
315 C = zeros(1,model.dim_x);
318 % Ccomp = {ones(1,model.dim_x)};
319 %Ccomp = {[2,0,1],[0,2,1]};
321 function Ccomp = C_components_expectation_M(model,model_data)
322 C = model_data.Ms(:)';
325 function Ccomp = C_components_expectation_N(model,model_data)
326 C = model_data.Ns(:)';
329 function Ccoeff = C_coefficients(model)
331 %Ccoeff = [model.C_weight,1-model.C_weight];
333 function D = D_func(model,model_data)
334 %D = [0;0]; % no decomposition required for scheme
335 D = [0]; % no decomposition required for scheme
337 function x0 = x0_func(model,model_data)
338 x0 = eval_affine_decomp(@x0_components,@x0_coefficients,model,model_data);
340 function x0comp = x0_components(model,model_data)
341 %if model.init_range == 0
342 % x0comp = {eye(model.dim_x,1)};
344 i = ((model_data.Ms+model_data.Ns)<=model.init_range);
345 %i = (model_data.Ms<=1);
346 x0comp = {i(:)/sum(i)}; % sum should be 1!
349 function x0coeff = x0_coefficients(model)
352 function model_data = my_ds_gen_model_data(model)
353 s = 1; % state counter
354 % generate list of states X and output weight vector
355 Ms = zeros(1,(model.range+1)*(model.range+1));
356 Ns = zeros(1,(model.range+1)*(model.range+1));
365 p = s-1; % number of states
367 %model_data = lin_ds_gen_model_data(model);
372 if model.affinely_decomposed
373 model.decomp_mode = 1;
374 % the following call is extremely expensive: 30 seconds!!!!
375 % therefore caching is activated
376 % model_data.A_components = model.A_function_ptr(model,model_data);
378 model.A_function_ptr,model,model_data);
379 model_data.B_components = model.B_function_ptr(model,model_data);
380 model_data.C_components = model.C_function_ptr(model,model_data);
382 % model_data.D_components = model.D_function_ptr(model);
383 model_data.x0_components = model.x0_function_ptr(model,model_data);
385 model_data.G = model.G_matrix_function_ptr(model,[]);
387 params.xnumintervals = model.range+1;
388 params.ynumintervals = model.range+1;
389 params.xrange = [-0.5,model.range+0.5];
390 params.yrange = [-0.5,model.range+0.5];
394 model_data.plot_grid =
rectgrid(params);
395 model_data.plot_nm_index_vector = ...
396 model_data.Ms + model_data.Ns*(model.range+1)+1;
398 function detailed_data = my_ds_gen_detailed_data(model,model_data)
399 detailed_data = lin_ds_gen_detailed_data(model,model_data);
401 detailed_data.plot_grid = model_data.plot_grid;
402 detailed_data.plot_nm_index_vector = model_data.plot_nm_index_vector;
404 function p = plot_state_probabilities(model,model_data,sim_data,params)
405 % plot of trajectory and output
406 %p = lin_ds_plot_sim_data(model,model_data,sim_data,params);
410 U = nan((model.range+1)*(model.range+1),size(sim_data.X,2));
411 U(model_data.plot_nm_index_vector,1:end)= sim_data.X;
413 %U(nm_index_vector,1) = model_data.Ns;
414 %U(nm_index_vector,2) = model_data.Ms;
415 %params.plot_function = @plot_element_data;
416 merged_plot_params = merge_model_plot_params(model, params);
417 merged_plot_params.plot_function = 'plot_element_data';
418 plot_data_sequence(model,model_data.plot_grid,U,merged_plot_params)
422 if ~isfield(model,'plot_title')
423 title('state probabilities');
425 title(model.plot_title);
428 function p = my_plot_slice(model,model_data,sim_data,params)
429 U = nan((model.range+1)*(model.range+1),1);
430 U(model_data.plot_nm_index_vector)= sim_data.X(:,params.timestep+1);
432 %U(nm_index_vector,1) = model_data.Ns;
433 %U(nm_index_vector,2) = model_data.Ms;
434 %params.plot_function = @plot_element_data;
435 params.plot_function = 'plot_element_data';
436 %colormap = jet(2^20);
437 params.colormap = my_transformed_colormap;
439 if ~isfield(params,'clim') && isfield(model,'clim')
440 params.clim = model.clim;
442 p = plot_element_data(model_data.plot_grid,U,params);
447 title('state probabilities');
449 function cm = my_transformed_colormap;
450 colormap = jet(65536);
452 %ind = round((
double(0:65535).^gamma)/(65535^gamma)*(2^20-1)+1);
453 ind = round((
double(0:255).^gamma)/(255^gamma)*65535+1);
454 cm = colormap(ind,:);
456 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
457 % CODE FROM STEFFEN, RESP HANWEI LI:
458 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
463 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
464 % $$$ % define parameter
471 % $$$ u1 = 0.01+(i-6)*0.001;
472 % $$$ u2 = 0.01+(i-6)*0.0005;
473 % $$$ k1 = 4*u1+(i-6)*0.0001;
476 % $$$ parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
480 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
481 % $$$ % build the nominal model
484 % $$$ s = 1; % state counter
486 % $$$ % generate list of states X
489 % $$$ for n=0:(range-m)
490 % $$$ X(:,s) = [m;n];
494 % $$$ p = s-1; % number of states
496 % $$$ % defining propensity functions
497 % $$$ a1 = @(x) (parameterlist(inomi,1)+parameterlist(inomi,4)*x(2)^h/(M1^h+x(2)^h));
498 % $$$ a2 = @(x) parameterlist(inomi,2)*x(1);
499 % $$$ a3 = @(x) parameterlist(inomi,5)*x(1)^h/(M2^h+x(1)^h);
500 % $$$ a4 = @(x) parameterlist(inomi,3)*x(2);
502 % $$$ A = sparse([],[],[],p,p,7*p);
504 % $$$ A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
505 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
506 % $$$ A(j,m) = a1(X(:,m));
507 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
508 % $$$ A(j,m) = a2(X(:,m));
509 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
510 % $$$ A(j,m) = a3(X(:,m));
511 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
512 % $$$ A(j,m) = a4(X(:,m));
516 % $$$ % KoeffizientenZerlegung
517 % $$$ clear A1; % for k1
518 % $$$ A1 = sparse([],[],[],p,p,7*p);
521 % $$$ j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
525 % $$$ clear A2; % for u1
526 % $$$ A2 = sparse([],[],[],p,p,7*p);
528 % $$$ A2(m,m) = -X(1,m);
529 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
530 % $$$ A2(j,m) = X(1,m);
533 % $$$ clear A3; % for u2
534 % $$$ A3 = sparse([],[],[],p,p,7*p);
536 % $$$ A3(m,m) = -X(2,m);
537 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
538 % $$$ A3(j,m) = X(2,m);
541 % $$$ clear A4; % for V1
542 % $$$ A4 = sparse([],[],[],p,p,7*p);
544 % $$$ A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
545 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
546 % $$$ A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
549 % $$$ clear A5; % for V2
550 % $$$ A5 = sparse([],[],[],p,p,7*p);
552 % $$$ A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
553 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
554 % $$$ A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
557 % $$$ % state space of the nominal system
561 % $$$ D = zeros(1,p);
562 % $$$ sys=ss(A,B,C,D);
564 % $$$ %Model reduction, sysb is a balanced realization for sys.sysb=T*sys*T^-1
565 % $$$ [sysb,g,T,Ti]= balreal(sys);
567 % $$$ sysr = modred(sysb,[n:p],'del');
569 % $$$ a2 = zeros(p-n+1,n-1);
570 % $$$ V = T^-1*[a1;a2];
571 % $$$ W = T'*[a1;a2];
573 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
574 % $$$ % compare the reduced system with different parameters with the original system
577 % $$$ % Ar = k1*A1r+u1*A2r+u2*A3r+V1*A4r+V2*A5r
578 % $$$ % Air = W'*Ai*V
579 % $$$ Ar = parameterlist(i,1)*W'*A1*V + parameterlist(i,2)*W'*A2*V + parameterlist(i,3)*W'*A3*V + parameterlist(i,4)*W'*A4*V + parameterlist(i,5)*W'*A5*V;
582 % $$$ % original system
583 % $$$ a1 = @(x) (parameterlist(i,1)+parameterlist(i,4)*x(2)^h/(M1^h+x(2)^h));
584 % $$$ a2 = @(x) parameterlist(i,2)*x(1);
585 % $$$ a3 = @(x) parameterlist(i,5)*x(1)^h/(M2^h+x(1)^h);
586 % $$$ a4 = @(x) parameterlist(i,3)*x(2);
588 % $$$ A = sparse([],[],[],p,p,7*p);
590 % $$$ A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
591 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
592 % $$$ A(j,m) = a1(X(:,m));
593 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
594 % $$$ A(j,m) = a2(X(:,m));
595 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
596 % $$$ A(j,m) = a3(X(:,m));
597 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
598 % $$$ A(j,m) = a4(X(:,m));
602 % $$$ subplot(4,3,i);
605 % $$$ [t,x]=ode15s(@(t,x) A*x, [0 1000000], x0);
606 % $$$ plot(t,sum(x'),'r');
610 % $$$ % reduced system
612 % $$$ x2=x1(1:(n-1),:);
613 % $$$ [tr,xr]=ode15s(@(tr,xr) Ar*xr, [0 1000000], x2);
616 % $$$ for m = 1:size(xr,1)
617 % $$$ y(m)=sysr.c*xr(m,:)';
621 % $$$ xlabel('zeit t');
function model = follicle_rect_model(params)
model of the human follicle growth
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
function [ RBext , dummy ] = RB_extension_PCA_fixspace(model, detailed_data)
function computing a RB basis extension for given parameters by the POD-Greedy algorithm.
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
function model = follicle_model(params)
model of the human follicle growth
A cartesian rectangular grid in two dimensions with axis parallel elements.
function save_detailed_simulations(model, model_data, M, savepath)
perform loop over detailed simulations and save results or check consistency with existing saved resu...