Version 0.6 generalized from two state unfolding to arbitrary
$N$-state rate reactions. This allows simulations of
-unfolding-refolding, metastable intermediates, etc.
+unfolding-refolding, metastable intermediates, etc., but required
+reorganizing the command line interface.
<<version definition>>=
#define VERSION "0.6"
@
\section{Main}
+\label{sec.main}
The unfolding system is stored as a chain of \emph{domains} (Figure
\ref{fig.domain_chain}). Domains can be folded, globular protein
The unfolding of the chain domains is modeled in two stages. First
the tension of the chain is calculated. Then the tension (assumed to
be constant over the length of the chain, see Section
-\ref{sec.timescales}), is applied to each folded domain, and
-used to calculate the probability of that subdomain unfolding in the
-next timestep $dt$. Then the time is stepped forward, and the process
-is repeated.
-
+\ref{sec.timescales}), is applied to each domain, and used to
+calculate the probability of that domain changing state in the next
+timestep $dt$. Then the time is stepped forward, and the process is
+repeated. The simulation is complete when a pre-set time limit
+[[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
<<main program>>=
int main(int argc, char **argv)
{
<<tension model getopt array definition>>
<<k model getopt array definition>>
- list_t *domain_list=NULL; /* variables initialized in get_options() */
+
+ /* variables initialized in get_options() */
+ list_t *state_list=NULL, *transition_list=NULL;
environment_t env={0};
- double P, dt_max, v, xstep;
- int num_folded=0, unfolding;
- double x=0, dt, F; /* variables used in the simulation loop */
- get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
- tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
+ double P, t_max, dt_max, v, x_step;
+ state_t *stop_state=NULL;
+
+ /* variables used in the simulation loop */
+ double t=0, x=0, dt, F; /* time, distance, timestep, force */
+ int transition=1; /* boolean marking a transition for tension optimization */
+
+ get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
+ NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
+ &state_list, &transition_list);
setup();
+
if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
- if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
- while (num_folded > 0) {
- F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
- if (xstep == 0)
- dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
+ if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
+ while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
+ F = find_tension(state_list, &env, x, transition);
+ if (x_step == 0)
+ dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
else
- dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0);
- unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
+ dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
+ transition=random_transitions(transition_list, F, dt, &env);
if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
- if (xstep == 0) {
+ t += dt;
+ if (x_step == 0) {
x += v*dt;
} else {
if (dt == dt_max) { /* step completed */
- x += xstep;
- dt_max = xstep / v;
+ x += x_step;
+ dt_max = x_step / v;
} else { /* still working on this step */
dt_max -= dt;
}
}
}
- destroy_domain_list(domain_list);
+ destroy_state_list(state_list);
+ destroy_transition_list(transition_list);
free_accels();
return 0;
}
@ The meat of the simulation is bundled into the three functions
-[[find_tension]], [[determine_dt]], and [[random_unfoldings]].
+[[find_tension]], [[determine_dt]], and [[random_transitions]].
[[find_tension]] is discussed in Section \ref{sec.find_tension},
[[determine_dt]] in Section \ref{sec.adaptive_dt}, and
-[[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
+[[random_transitions]] in Section \ref{sec.transition_rate}.
The stretched distance is extended in one of two modes depending on
[[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
calculated, bridging the gap between theory and experiment.
Environmental parameters important in determining reaction rates and
-tensions (e.g. temperature, pH) are stored in a single structure to
+tensions (e.g.\ temperature, pH) are stored in a single structure to
facilitate extension to more complicated models in the future.
<<environment definition>>=
typedef struct environment_struct {
<<find tension>>
<<determine dt>>
<<happens>>
-<<does domain unfold>>
-<<randomly unfold domains>>
+<<does domain transition>>
+<<randomly transition domains>>
@
\subsection{Tension}
Because the stretched system may be made up of several parts (folded
domains, unfolded domains, spring-like cantilever, \ldots), we will
-assign the domains to models and groups. The models are used to
-determine a domain's tension (Hookean spring, WLC, \ldots). Groups
-will represent collections of domains which the model should treat as
-a single entity. A domain may belong to different groups or models
-depending on its state. For example, a domain might be modeled as a
-freely-jointed chain when it is folded and as a worm-like chain when
-it is unfolded. The domains are assumed to be commutative, so
-ordering is ignored. The interactions between the groups are assumed
-to be linear, but the interactions between domains of the same group
-need not be. This allows for non-linear group models such as th
-worm-like or freely-jointed chains. Each model has a tension handler
-function, which gives the tension $F$ needed to stretch a given group
-of domains a distance $x$.
-
-To understand the purpose of groups, consider a chain of proteins
-where two folded proteins $A$ and $B$ are modeled as WLCs with
-persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
-modeled as WLCs with persistence length $p_u$. The proteins are
-attached to a cantilever $E$ qof spring constant $κ$. The string
-would get split into three lists:
-\begin{center}
- \begin{tabular}{llll}
- Model & Group & List & Tension \\
- WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
- WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
- Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
- \end{tabular}
-\end{center}
-Note that group numbers only matter \emph{within} model classes, since
-grouping domains with seperate models doesn't make sense.
+assign the domains to states with tension models and groups. The
+models are used to determine the tension of all the domain in a given
+state following some model (Hookean spring, WLC, \ldots). The domains
+are assumed to be commutative, so ordering is ignored. The
+interactions between the groups are assumed to be linear, but the
+interactions between domains of the same group need not be. Each
+model has a tension handler function, which gives the tension $F$
+needed to stretch a given group of domains a distance $x$.
+
+GROUPS ARE CURRENTLY DISABLED.
+
+Groups represent collections of states which the model should treat as
+a single entity. To understand the purpose of groups, consider a
+system with two unfolded domain states (e.g.\ for two different
+proteins) that were both modeled as WLCs. If both unfolded states had
+the same persistence length, it would be wise to assign them to the
+same group. This would reduce both the computational cost of
+balancing the tension and the error associated with the inter-group
+interaction linearity. Note that group numbers only matter
+\emph{within} model classes, since grouping domains with seperate
+models doesn't make sense.
<<find tension definitions>>=
#define NUM_TENSION_MODELS 5
<<tension handler types>>=
typedef struct tension_handler_data_struct {
- /* int num_domains; */
- /* domain_t *domains; */
- list_t *group;
+ list_t *group_tension_params; /* one item for each domain in the group */
environment_t *env;
void *persist;
} tension_handler_data_t;
After sorting the chain into separate groups $G_i$, with tension
handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
-\forall i,j$. Note that there may be several groups within each model
-class. [[find_tension]] is basically a wrapper to feed proper input
-into the [[tension_balance]] function. [[unfolding]] should be set to
-0 if there were no unfoldings since the previous call to
-[[find_tension]]. While were messing with the tension handlers, we
-also grab a measure of the current linker stiffness for the
-environmental variable, which is needed by some of the unfolding rate
-models (e.g. the linker-stiffness corrected Bell model, Appendix
+\forall i,j$. Note that there may be several states within each
+group. [[find_tension]] is basically a wrapper to feed proper input
+into the [[tension_balance]] function. [[transition]] should be set
+to 0 if there were no transitions since the previous call to
+[[find_tension]] to support some optimizations that assume a small
+increase in tension between steps (only true if no transition has
+occured). While were messing with the tension handlers, we also grab
+a measure of the current linker stiffness for the environmental
+variable, which is needed by some of the unfolding rate models
+(e.g.\ the linker-stiffness corrected Bell model, Appendix
\ref{sec.kbell}).
<<find tension>>=
-double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
+double find_tension(list_t *state_list, environment_t *env, double x, int transition)
{
static int num_active_groups=0;
static one_dim_fn_t **per_group_handlers = NULL;
- static void **per_group_params = NULL;
+ static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
static double last_x;
tension_handler_data_t *tdata;
double F;
int i;
#ifdef DEBUG
- fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
- list_print(stderr, domain_list, "domain list");
+ fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
#endif
- if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
- /* free old statics */
- if (per_group_handlers != NULL)
- free(per_group_handlers);
- if (per_group_params != NULL) {
- for (i=0; i < num_active_groups; i++) {
- assert(per_group_params[i] != NULL);
- free(per_group_params[i]);
- }
- free(per_group_params);
- }
-
- /* get new statics */
- get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
+ if (transition != 0 || num_active_groups == 0) { /* setup statics */
+ /* get new statics, freeing the old ones if they aren't NULL */
+ get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_data);
/* fill in the group handlers and params */
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
+#endif
for (i=0; i<num_active_groups; i++) {
- tdata = (tension_handler_data_t *) per_group_params[i];
+ tdata = (tension_handler_data_t *) per_group_data[i];
#ifdef DEBUG
- fprintf(stderr, "active group %d of %d", i, num_active_groups);
- list_print(stderr, tdata->group, " parameters");
+ fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
+ list_print(stderr, tdata->group_tension_params, " parameters");
#endif
tdata->env = env;
/* tdata->persist continues unchanged */
}
last_x = -1.0;
- } /* else roll with the current statics */
+ } /* else continue with the current statics */
- F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
-
- {
- double dx;
- double xlow;
- double Flow;
-#ifdef DEBUG
- fprintf(stderr, "finding linker stiffness at distance %g, tension %g\n", x, F);
-#endif
- if (last_x == -1.0)
- dx = x/200.0;
- else
- dx = (x-last_x);
- if (dx == 0) { /* e.g. if x == 0 */
- env->stiffness = 0;
- } else {
- xlow = x-dx;
- Flow = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, xlow);
- env->stiffness = (F-Flow)/dx;
- }
-#ifdef DEBUG
- fprintf(stderr, " stiffness (%g-%g)/%g = %g\t(dx = %g = %g-%g)\n", F, Flow, dx, env->stiffness, dx, x, xlow);
-#endif
- }
+ F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
last_x = x;
ourselves to strictly monotonically increasing, non-negative tension
functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
simpler. For example, it guarantees the existence of a unique, real
-solution for finite forces. See Appendix \ref{app.tension_balance}
+solution for finite forces. See Appendix \ref{sec.tension_balance}
for the tension balancing implementation.
Each group must define a parameter structure if the tension model is
a function to destroy the parameter structure, and
a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
The tension-specific parameter structures are contained in the domain
-group field. See the Section \ref{app.model_selection} for a
+group field. See the Section \ref{sec.model_selection} for a
discussion on the command line framework. See the worm-like chain
implementation for an example (Section \ref{sec.wlc}).
are Markovian (memory-less), which is not really a problem for the
chains (see \citet{evans99} for WLC thermalization rate calculations).
-\subsection{Unfolding rate}
-\label{sec.unfolding_rate}
+\subsection{Transition rate}
+\label{sec.transition_rate}
-Each folded domain is modeled as unimolecular, first order reaction
+Each state transition is modeled as unimolecular, first order reaction
$$
- \text{Folded} \xrightarrow{k} % in tex: X atop Y
- \text{Unfolded}
+ \text{State 1} \xrightarrow{k} \text{State 2}
$$
With the rate constant $k$ defined by
$$
- \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
+ \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
$$
-So the probability of a given protein unfolding in a short time $dt$
-is given by
+So the probability of a given domain transitioning in a short time
+$dt$ is given by
$$
P = k \cdot dt.
$$
-<<does domain unfold>>=
-int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
-{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
+<<does domain transition>>=
+int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
+{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
double k;
- k = accel_k(domain->k, F, env, domain->k_params);
- //(*domain->k)(F, env, domain->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
+ //(*transition->k)(F, env, domain->k_params);
//printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
- return happens(k*dt); /* dice roll for prob. k*dt event */
-}
-@ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
-
-<<randomly unfold domains>>=
-int random_unfoldings(list_t *domain_list, double tension,
- double dt, environment_t *env,
- int *num_folded_domains)
-{ /* returns 1 if there was an unfolding and 0 otherwise */
- while (domain_list != NULL) {
- if (D(domain_list)->state == DS_FOLDED
- && domain_unfolds(tension, dt, env, D(domain_list))) {
- if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
- fprintf(stdout, "%g\n", tension);
- D(domain_list)->state = DS_UNFOLDED;
- (*num_folded_domains)--;
- return 1; /* our one domain has unfolded, stop looking for others */
+ return happens(k*dt*num_domains); /* N dice rolls for prob. k*dt event */
+}
+@ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
+
+<<randomly transition domains>>=
+int random_transitions(list_t *transition_list, double tension,
+ double dt, environment_t *env)
+{ /* returns 1 if there was a transition and 0 otherwise */
+ while (transition_list != NULL) {
+ if (T(transition_list)->initial_state->num_domains > 0
+ && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
+ if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
+ fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
+ T(transition_list)->initial_state->num_domains--;
+ T(transition_list)->final_state->num_domains++;
+ return 1; /* our one domain has transitioned, stop looking for others */
}
- domain_list = domain_list->next;
+ transition_list = transition_list->next;
}
return 0;
}
\label{sec.adaptive_dt}
We'd like to pick $dt$ so the probability of changing state
-(e.g. unfolding) in the next timestep is small. If we don't adapt our
+(e.g.\ unfolding) in the next timestep is small. If we don't adapt our
timestep, we also risk breaking our approximation that $F(x' \in [x,
x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
given timestep. The simulation would have been accurate for
The actual adaptive timestep implementation is not particularly
interesting, since we are only required to reduce $dt$ to somewhere
below a set threshold, so I've removed it to Appendix
-\ref{app.adaptive_dt}.
+\ref{sec.adaptive_dt_implementation}.
\section{Models}
<<help>>
<<index tension model>>
<<index k model>>
-<<generate domain>>
<<get options>>
@
\subsection{Model selection}
-\label{app.model_selection}
+\label{sec.model_selection}
The main difficulty with the command line interface in \prog\ is
developing an intuitive method for accessing the possibly large number
<<k model getopt definitions>>
@
+In order to choose models, we need a method of assembling a reaction
+graph such as that in Figure \ref{fig.domain_states} and population
+the graph with a set of domains. First we declare the domain states
+following
+\begin{center}
+[[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
+\end{center}
+e.g.
+\begin{center}
+[[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
+\end{center}
+This creates the state named [[unfolded]], using the WLC model and one
+as a group number (defaults to 0, see Section \ref{sec.find_tension}).
+The tension parameters are then set to [[1e-8,4e-10]], which in the
+case of the WLC are the contour and persistence lengths respectively.
+Finally we populate the state by adding five domains to the state.
+The name of the state is importent for identifying states later.
+
+Once the domains have been created and populated, you can added
+transition rates following
+\begin{center}
+ [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
+\end{center}
+e.g.
+\begin{center}
+ [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
+\end{center}
+This creates a reaction going from the [[folded]] state to the
+[[unfolded]] state, with the rate constant given by the Bell model
+(Appendix \ref{sec.bell}). The unfolding parameters are then set to
+[[3.3e-4,0.25e-9]], which in the case of the Bell model are the
+unforced rate and transition state distance respectively.
\subsubsection{Tension}
model.
<<tension model getopt definitions>>=
typedef struct tension_model_getopt_struct {
- one_dim_fn_t *handler; /* fn implementing the model on a group */
+ one_dim_fn_t *handler; /* fn implementing the model on a group */
char *name; /* name string identifying the model */
char *description; /* a brief description of the model */
int num_params; /* number of extra params for the model */
\subsection{help}
<<help>>=
-void help(char *prog_name, double P, double dt_max, double v, double xstep,
- environment_t *env,
+void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
+ state_t *stop_state, environment_t *env,
int n_tension_models, tension_model_getopt_t *tension_models,
- int folded_tension_model, int unfolded_tension_model,
int n_k_models, k_model_getopt_t *k_models,
- int k_model)
+ int k_model, list_t *state_list)
{
+ state_t *state=NULL;
int i, j;
+
+ if (state_list != NULL)
+ state = S(tail(state_list));
+
printf("usage: %s [options]\n", prog_name);
printf("Version %s\n\n", VERSION);
printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
+
printf("Simulation options:\n");
+ printf("");
+ printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
+ printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
+ printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
printf("-P P\tTarget probability for dt (currently %g)\n", P);
- printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
- printf("-x dx\tPulling step size (currently %g m)\n", xstep);
- printf("\tWhen dx = 0, the pulling is continuous\n");
- printf("\tWhen dx > 0, the pulling occurs in discrete steps\n");
+ printf("-x dx\tPulling step size (currently %g m)\n", x_step);
+ printf("\tWhen dx = 0, the pulling is continuous.\n");
+ printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
+
printf("Environmental options:\n");
printf("-T T\tTemperature T (currently %g K)\n", env->T);
printf("-C T\tYou can also set the temperature T in Celsius\n");
- printf("Model options:\n");
- printf("The domains exist in either the folded or unfolded state\n");
- printf("The following options change the default behavior in each state.\n");
- printf("-m model[,group]\tFolded tension model (currently %s)\n",
- tension_models[folded_tension_model].name);
- printf("-a args\tFolded tension model argument string (currently %s)\n",
- tension_models[folded_tension_model].params);
- printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
- tension_models[unfolded_tension_model].name);
- printf("-A args\tUnfolded tension model argument string (currently %s)\n",
- tension_models[unfolded_tension_model].params);
- printf("The following options change the unfolding rate.\n");
- printf("-k model\tTransition rate model (currently %s)\n",
- k_models[k_model].name);
- printf("-K args\tTransition rate model argument string (currently %s)\n",
- k_models[k_model].params);
- printf("Domain creation options:\n");
- printf("Once you've set up the appropriate default models, you need to add the domains\n");
- printf("-F n\tAdd n folded domains with the current models\n");
- printf("-U n\tAdd n unfolded domains with the current models\n");
- printf("Output mode options:\n");
- printf("There are two output modes. In standard mode, only the unfolding\n");
+
+ printf("State creation options:\n");
+ printf("-s name,model[[,group],params]\tCreate a new state.\n");
+ printf("Once you've set up the state, you may populate it with domains\n");
+ printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
+
+ printf("Transition creation options:\n");
+ printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
+
+ printf("To double check your reaction graph, you can produce graphviz dot output\n");
+ printf("describing the current states and transitions.\n");
+ printf("-D\tPrint dot description of states and transitions and exit.\n");
+
+ printf("Simulation output mode options:\n");
+ printf("There are two output modes. In standard mode, only the transition\n");
printf("events are printed. For example:\n");
- printf(" #Unfolding Force (N)\n");
- printf(" 123.456e-12\n");
+ printf(" #Force (N)\tInitial state\tFinal state\n");
+ printf(" 123.456e-12\tfolded\tunfolded\n");
printf(" ...\n");
printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
- printf(" 0.001\t0.0023\n");
+ printf(" 0.001\t0.0023\t0.002\n");
printf(" ...\n");
- printf("-V\tChange output to verbose mode\n");
- printf("-h\tPrint this help and exit\n");
+ printf("-V\tChange output to verbose mode.\n");
+ printf("-h\tPrint this help and exit.\n");
printf("\n");
+
printf("Tension models:\n");
for (i=0; i<n_tension_models; i++) {
printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
printf("\t\tdefaults: %s\n", k_models[i].params);
}
+
printf("Examples:\n");
- printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
- printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
- printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
- printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
+ printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
+ printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s 'unfolded,wlc,{0.39e-9,28e-9}' -k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' -q folded\n", prog_name);
+ printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
+ printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s 'folded,wlc,{1e-8,4e-10}' -N8 -s 'unfolded,wlc,1,{0.4e-9,28.1e-9}' -k 'folded,unfolded,bell,{5e-5,0.225e-9}' -t0.25\n", prog_name);
}
@
<<get options>>=
void get_options(int argc, char **argv,
- double *pP, double *pDt_max, double *pV, double *pXstep,
- environment_t *env,
+ double *pP, double *pT_max, double *pDt_max, double *pV,
+ double *pX_step, state_t **pStop_state, environment_t *env,
int n_tension_models, tension_model_getopt_t *tension_models,
int n_k_models, k_model_getopt_t *k_models,
- list_t **pDomain_list, int *pNum_folded)
+ list_t **pState_list, list_t **pTransition_list)
{
char *prog_name = NULL;
- char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh";
- int ftension_model=0, utension_model=0, k_model=0;
- char *ftension_params=NULL, *utension_params=NULL;
- int i, n, ftension_group=0, utension_group=0;
+ char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
+ int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
+ char *state_name=NULL;
+ int i, n, tension_group=0;
+ list_t *temp_list=NULL;
int num_strings;
char **strings;
+ void *model_params; /* for INIT_MODEL() */
+ int num_param_args; /* for INIT_MODEL() */
+ char **param_args; /* for INIT_MODEL() */
extern char *optarg;
extern int optind, optopt, opterr;
#endif
/* setup defaults */
- flags = FLAG_OUTPUT_UNFOLDING_FORCES;
+ flags = FLAG_OUTPUT_TRANSITION_FORCES;
prog_name = argv[0];
- *pP = 1e-3; /* % pop per s */
- *pDt_max = 0.001; /* s */
- *pV = 1e-6; /* m/s */
- *pXstep = 0.0; /* m */
- env->T = 300.0; /* K */
+ *pStop_state = NULL;
+ *pT_max = -1; /* s, -1 removes this limit */
+ *pDt_max = 0.001; /* s */
+ *pP = 1e-3; /* % pop per s */
+ *pV = 1e-6; /* m/s */
+ *pX_step = 0.0; /* m */
+ env->T = 300.0; /* K */
+ *pState_list = NULL;
+ *pTransition_list = NULL;
while ((c=getopt(argc, argv, options)) != -1) {
switch(c) {
+ case 'q':
+ if (STRMATCH(optarg, "-")) {
+ *pStop_state = NULL;
+ } else {
+ temp_list = head(*pState_list);
+ while (temp_list != NULL) {
+ if (STRMATCH(S(temp_list)->name, optarg)) {
+ *pStop_state = S(temp_list);
+ break;
+ }
+ temp_list = temp_list->next;
+ }
+ }
+ break;
+ case 't': *pT_max = atof(optarg); break;
+ case 'd': *pDt_max = atof(optarg); break;
case 'P': *pP = atof(optarg); break;
- case 't': *pDt_max = atof(optarg); break;
case 'v': *pV = atof(optarg); break;
- case 'x': *pXstep = atof(optarg); break;
+ case 'x': *pX_step = atof(optarg); break;
case 'T': env->T = atof(optarg); break;
case 'C': env->T = atof(optarg)+273.15; break;
- case 'm':
+ case 's':
parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
- assert(num_strings == 1 || num_strings == 2);
- if (num_strings == 1)
- ftension_group = 0;
+ assert(num_strings >= 2 && num_strings <= 4);
+
+ state_name = strings[0];
+ if (state_index(*pState_list, state_name) != -1) {
+ fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
+ exit(1);
+ }
+ tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
+ if (num_strings == 4)
+ tension_group = atoi(strings[2]);
else
- ftension_group = atoi(strings[1]);
- ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
- ftension_params = tension_models[ftension_model].params;
+ tension_group = 0;
+ if (num_strings >= 3)
+ tension_models[tension_model].params = strings[num_strings-1];
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s (%s), group %g\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
+#endif
+ INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
+ push(pState_list, create_state(state_name,
+ tension_models[tension_model].handler,
+ tension_group,
+ model_params,
+ tension_models[tension_model].destructor,
+ 0), 1);
+ *pState_list = tail(*pState_list);
+ state_name = NULL;
free_parsed_list(num_strings, strings);
break;
- case 'a':
- ftension_params = optarg;
+ case 'N':
+ n = atoi(optarg);
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
+#endif
+ S(*pState_list)->num_domains += n;
break;
- case 'M':
+ case 'k':
parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
- assert(num_strings == 1 || num_strings == 2);
- if (num_strings == 1)
- utension_group = 0;
- else
- utension_group = atoi(strings[1]);
- utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
- utension_params = tension_models[utension_model].params;
+ assert(num_strings >= 3 && num_strings <= 4);
+ initial_state = state_index(*pState_list, strings[0]);
+ final_state = state_index(*pState_list, strings[1]);
+ k_model = index_k_model(n_k_models, k_models, strings[2]);
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t%s -> %s reaction from state %s (%d) -> state %s (%d)\n", __FUNCTION__, optarg, strings[2], strings[0], initial_state, strings[1], final_state);
+#endif
+ assert(initial_state != final_state);
+ if (num_strings == 4)
+ k_models[k_model].params = strings[3];
+ INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
+ push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
+ list_element(*pState_list, final_state),
+ k_models[k_model].k,
+ model_params,
+ k_models[k_model].destructor), 1);
free_parsed_list(num_strings, strings);
break;
- case 'A':
- utension_params = optarg;
- break;
- case 'k':
- k_model = index_k_model(n_k_models, k_models, optarg);
- break;
- case 'K':
- k_models[k_model].params = optarg;
- break;
- case 'F':
- n = atoi(optarg);
- assert(n > 0);
- for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_FOLDED,
- tension_models+ftension_model,
- ftension_group,
- ftension_params,
- tension_models+utension_model,
- utension_group,
- utension_params,
- k_models+k_model), 1);
- }
- *pNum_folded += n;
- break;
- case 'U':
- n = atoi(optarg);
- assert(n > 0);
- for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_UNFOLDED,
- tension_models+ftension_model,
- ftension_group,
- ftension_params,
- tension_models+utension_model,
- utension_group,
- utension_params,
- k_models+k_model), 1);
- }
+ case 'D':
+ print_state_graph(stdout, *pState_list, *pTransition_list);
+ exit(0);
break;
case 'V':
flags = FLAG_OUTPUT_FULL_CURVE;
fprintf(stderr, "unrecognized option '%c'\n", optopt);
/* fall through to default case */
default:
- help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
+ help(prog_name, *pP, *pT_max, *pDt_max, *pV, *pX_step, *pStop_state, env, n_tension_models, tension_models, n_k_models, k_models, k_model, *pState_list);
exit(1);
}
}
/* check the arguments */
- assert(*pDomain_list != NULL); /* no domains! */
+ assert(*pState_list != NULL); /* no domains! */
assert(*pP > 0.0); assert(*pP < 1.0);
assert(*pDt_max > 0.0);
assert(*pV > 0.0);
assert(env->T > 0.0);
+
+ crosslink_groups(*pState_list, *pTransition_list);
+
return;
}
@
*/
#define INIT_MODEL(role, model, param_string, param_pointer) \
do { \
- parse_list_string(param_string, SEP, '{', '}', \
+ parse_list_string((param_string), SEP, '{', '}', \
&num_param_args, ¶m_args); \
- if (num_param_args != model->num_params) { \
+ if (num_param_args != (model)->num_params) { \
fprintf(stderr, \
"%s model %s expected %d params,\n", \
- role, model->name, model->num_params); \
+ (role), (model)->name, (model)->num_params); \
fprintf(stderr, \
"not the %d params in '%s'\n", \
- num_param_args, param_string); \
- assert(num_param_args == model->num_params); \
+ num_param_args, (param_string)); \
+ assert(num_param_args == (model)->num_params); \
} \
- if (model->creator) \
- param_pointer = (*model->creator)(param_args); \
+ if ((model)->creator) \
+ param_pointer = (*(model)->creator)(param_args); \
else \
param_pointer = NULL; \
free_parsed_list(num_param_args, param_args); \
} while (0);
@
-<<generate domain>>=
-void *generate_domain(enum domain_state_t state,
- tension_model_getopt_t *folded_model,
- int folded_group,
- char *folded_param_string,
- tension_model_getopt_t *unfolded_model,
- int unfolded_group,
- char *unfolded_param_string,
- k_model_getopt_t *k_model)
-{
- void *folded_params, *unfolded_params, *k_params;
- int num_param_args; /* for INIT_MODEL() */
- char **param_args; /* for INIT_MODEL() */
-
-#ifdef DEBUG
- fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
- fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
- k_model->name, k_model->params,
- folded_model->name, folded_model->handler, folded_group, folded_param_string,
- unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
-#endif
- INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
- INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
- INIT_MODEL("k", k_model, k_model->params, k_params);
- return create_domain(state,
- k_model->k, k_params, k_model->destructor,
- folded_model->handler, folded_group, folded_params, folded_model->destructor,
- unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
- );
-}
-@
-
\phantomsection
\appendix
\addcontentsline{toc}{section}{Appendicies}
The general layout of our simulation code is:
<<*>>=
+//#define DEBUG
<<license comment>>
<<includes>>
<<definitions>>
<<string comparison definition and globals>>
<<initialize model definition>>
<<get options definitions>>
-<<domain definitions>>
+<<state definitions>>
+<<transition definitions>>
@
<<globals>>=
@
<<functions>>=
-<<create/destroy domain>>
-<<destroy domain list>>
+<<create/destroy state>>
+<<destroy state list>>
+<<state index>>
+<<create/destroy transition>>
+<<destroy transition list>>
+<<print state graph>>
<<group functions>>
<<simulation functions>>
<<boilerplate functions>>
<<flag definitions>>=
/* in octal b/c of prefixed '0' */
-#define FLAG_OUTPUT_FULL_CURVE 01
-#define FLAG_OUTPUT_UNFOLDING_FORCES 02
+#define FLAG_OUTPUT_FULL_CURVE 01
+#define FLAG_OUTPUT_TRANSITION_FORCES 02
@
<<flag globals>>=
@
\subsection{Utilities}
-\label{app.utils}
+\label{sec.utils}
<<max/min definitions>>=
#define MAX(a,b) ((a)>(b) ? (a) : (b))
\subsection{Adaptive timesteps}
-\label{app.adaptive_dt}
-
-$F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
-so basing the timestep on the the unfolding probability at the current tension
-is dangerous, and we need to search for a $dt$ for which
-$P(F(x+v*dt)) < P_\text{target}$.
-There are two cases to consider.
-In the most common, no domains have unfolded since the last step,
-and we expect the next step to be slightly shorter than the previous one.
-In the less common, domains did unfold in the last step,
-and we expect the next step to be considerably longer than the previous one.
+\label{sec.adaptive_dt_implementation}
+
+$F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
+chain model), so basing the timestep on the the unfolding probability
+at the current tension is dangerous, and we need to search for a $dt$
+for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
+consider. In the most common, no domains have unfolded since the last
+step, and we expect the next step to be slightly shorter than the
+previous one. In the less common, domains did unfold in the last
+step, and we expect the next step to be considerably longer than the
+previous one.
<<search dt>>=
-double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
- list_t *domain_list,
+double search_dt(transition_t *transition,
+ list_t *state_list,
environment_t *env, double x,
double target_prob, double max_dt, double v)
-{ /* Returns the timestep dt in seconds for the current folded domain.
- * Takes a list of tension handlers, the list of domains,
- * a pointer env to the environmental data, a starting separation x in m,
- * a target_prob between 0 and 1,
+{ /* Returns a valid timestep dt in seconds for the given transition.
+ * Takes a list of all states, a pointer env to the environmental data,
+ * a starting separation x in m, a target_prob between 0 and 1,
* max_dt in s, stretching velocity v in m/s.
*/
double F, k, dtCur, dtU, dtUCur, dtL, dt;
/* get upper bound using the current position */
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
+ F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
//printf("Start with x = %g (F = %g)\n", x, F);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
//printf("x %g\tF %g\tk %g\n", x, F, k);
dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
if (dtU > max_dt) {
//printf("overshot max_dt\n");
dtU = max_dt;
}
- if (v == 0) /* discrete stepping, so force is temporatily constant */
+ if (v == 0) /* discrete stepping, so force is temporarily constant */
return dtU;
/* set a lower bound on dt too */
/* The dt determined above may produce illegitimate forces or ks.
* Reduce the upper bound until we have valid ks. */
dt = dtU;
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0);
while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
dtU /= 2.0;
dt = dtU;
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0);
}
//printf("Try for dt = %g (F = %g)\n", dt, F);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
/* returned k may be -1.0 */
//printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
while (k == -1.0) { /* reduce step until we hit a valid k */
dtU /= 2.0;
dt = dtU; /* hopefully, we can use the max dt, see if it works */
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0);
//printf("Try for dt = %g (F = %g)\n", dt, F);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
//printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
}
assert(dtU > 1e-14); /* timestep to valid k too small */
/* dtU wasn't safe, lets see what would be. */
while (dtU > 1.1*dtL) { /* until the range is reasonably small */
dt = (dtU + dtL) / 2.0;
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0);
//printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
dtCur = target_prob / k;
//printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
if (dtCur > dt) /* safe timestep back from x+dt covers dt */
}
return MAX(dtUCur, dtL);
}
+
@
To determine $dt$ for an array of potentially different folded domains,
we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
<<determine dt>>=
<<search dt>>
-double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
+double determine_dt(list_t *state_list, list_t *transition_list,
environment_t *env, double x,
double target_prob, double dt_max, double v)
{ /* Returns the timestep dt in seconds.
- * Takes the list of folded domains, target_prob between 0 and 1,
- * F in N, and T in K. */
+ * Takes the state and transition lists, a pointer to the environment,
+ * an extension x in m, target_prob between 0 and 1, dt_max in s,
+ * and v in m/s. */
double dt=dt_max, new_dt;
assert(target_prob > 0.0); assert(target_prob < 1.0);
assert(dt_max > 0.0);
-
- while (domain_list != NULL) {
- if (D(domain_list)->state == DS_FOLDED) {
- new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
+ assert(state_list != NULL);
+ assert(transition_list != NULL);
+ transition_list = head(transition_list);
+ while (transition_list != NULL) {
+ if (T(transition_list)->initial_state->num_domains > 0) {
+ new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
dt = MIN(dt, new_dt);
}
- domain_list = domain_list->next;
+ transition_list = transition_list->next;
}
return dt;
}
-@
-
-\subsection{Domain data}
-Currently domains exist in two states, folded and unfolded, and the
-only allowed transitions are folded $\rightarrow$ unfolded. Of
-course, it wouldn't be too complicated to extent this to a multi-state
-system, with an array containing the domains group for each possible
-state, and a matrix of transition-rate-calculating functions.
-However, at this point such generality seems unnecessary at this
-point.
-<<domain definitions>>=
-enum domain_state_t {DS_FOLDED,
- DS_UNFOLDED
-};
+@
-typedef struct domain_struct {
- enum domain_state_t state;
- one_dim_fn_t *folded_handler;
- int folded_group;
- one_dim_fn_t *unfolded_handler;
- int unfolded_group;
- k_func_t *k; /* function returning unfolding rate */
- void *folded_params; /* pointer to folded parameters */
- void *unfolded_params; /* pointer to unfolded parameters */
- void *k_params; /* pointer to k parameters */
- destroy_data_func_t *destroy_folded;
- destroy_data_func_t *destroy_unfolded;
- destroy_data_func_t *destroy_k;
-} domain_t;
-
-/* get the domain data for the current list node */
-#define D(list) ((domain_t *)(list)->d)
-/* get the tension params for the current list node */
-#define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
- ? ((domain_t *)(list)->d)->folded_params \
- : ((domain_t *)(list)->d)->unfolded_params)
-@
-[[k]] is a pointer to the function determining the unfolding rate for a given tension.
-[[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
-[[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
-The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
-We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
-
-[[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
-<<create/destroy domain>>=
-domain_t *create_domain(enum domain_state_t state,
- k_func_t *k,
- void *k_params,
- destroy_data_func_t *destroy_k,
- one_dim_fn_t *folded_handler,
- int folded_group,
- void *folded_params,
- destroy_data_func_t *destroy_folded,
- one_dim_fn_t *unfolded_handler,
- int unfolded_group,
- void *unfolded_params,
- destroy_data_func_t *destroy_unfolded)
-{
- domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
+\subsection{State data}
+\label{sec.state_data}
+
+The domains exist in one of an array of possible states. Each state
+has a tension model which determines it's force/extension curve
+(possibly [[null]] for rigid states, see Appendix
+\ref{sec.null_tension}).
+
+Groups are, for example, for two domain states with WLC tensions; one
+with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
+$L=28\U{nm}$. Since the persistence lengths are the same, you would
+like to add the contour lengths of all the domains in \emph{both}
+states, and plug that total contour length into a single WLC
+evaluation (see Section \ref{sec.find_tension}).
+<<state definitions>>=
+typedef struct state_struct {
+ char *name; /* name string */
+ one_dim_fn_t *tension_handler; /* tension handler */
+ int tension_group; /* for grouping similar tension states */
+ void *tension_params; /* pointer to tension parameters */
+ destroy_data_func_t *destroy_tension;
+ int num_domains; /* number of domains currently in state */
+ /* cached lookups generated from the above data */
+ int num_out_trans; /* length of out_trans array */
+ struct transition_struct **out_trans; /* transitions leaving this state */
+ int num_grouped_states; /* length of grouped states array */
+ struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
+} state_t;
+
+/* get the state data for the current list node */
+#define S(list) ((state_t *)(list)->d)
+
+@ [[tension_params]] is a pointer to the parameters used by the
+handler function when determining the tension. [[destroy_tension]]
+points to a function for freeing [[tension_params]]. We it in the
+state structure so that [[destroy_state]] and will have an easier job
+cleaning up.
+
+[[create_]] and [[destroy_state]] are simple wrappers around
+[[malloc]] and [[free]].
+<<create/destroy state>>=
+state_t *create_state(char *name,
+ one_dim_fn_t *tension_handler,
+ int tension_group,
+ void *tension_params,
+ destroy_data_func_t *destroy_tension,
+ int num_domains)
+{
+ state_t *ret = (state_t *)malloc(sizeof(state_t));
assert(ret != NULL);
- if (state == DS_FOLDED) {
- assert(k != NULL); /* the pointer points somewhere valid */
- assert(*k != NULL); /* and there is something useful there */
- } else
- assert(state == DS_UNFOLDED);
- ret->state = state;
- ret->folded_handler = folded_handler;
- ret->folded_group = folded_group;
- ret->unfolded_handler = unfolded_handler;
- ret->unfolded_group = unfolded_group;
- ret->k = k;
- ret->folded_params = folded_params;
- ret->unfolded_params = unfolded_params;
- ret->k_params = k_params;
- ret->destroy_folded = destroy_folded;
- ret->destroy_unfolded = destroy_unfolded;
- ret->destroy_k = destroy_k;
+ /* make a copy of the name, so we aren't dependent on the original */
+ ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
+ assert(ret->name != NULL);
+ strcpy(ret->name, name); /* we know ret->name is long enough */
+
+ ret->tension_handler = tension_handler;
+ ret->tension_group = tension_group;
+ ret->tension_params = tension_params;
+ ret->destroy_tension = destroy_tension;
+ ret->num_domains = num_domains;
+
+ ret->num_out_trans = 0;
+ ret->out_trans = NULL;
+ ret->num_grouped_states = 0;
+ ret->grouped_states = NULL;
+
#ifdef DEBUG
- fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
+ fprintf(stderr, "generated state %s (%p) with tension handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->tension_params, ret->tension_group);
#endif
return ret;
}
-void destroy_domain(domain_t *domain)
+void destroy_state(state_t *state)
{
- if (domain) {
- //printf("domain %p & %p\n", *domain, domain);
- if (domain->destroy_folded)
- (*domain->destroy_folded)(domain->folded_params);
- if (domain->destroy_unfolded)
- (*domain->destroy_unfolded)(domain->unfolded_params);
- if (domain->destroy_k)
- (*domain->destroy_k)(domain->k_params);
- free(domain);
+ if (state) {
+#ifdef DEBUG
+ fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
+#endif
+ if (state->name)
+ free(state->name);
+ if (state->destroy_tension)
+ (*state->destroy_tension)(state->tension_params);
+ if (state->out_trans)
+ free(state->out_trans);
+ if (state->grouped_states)
+ free(state->grouped_states);
+ free(state);
}
}
+
@
-<<destroy domain list>>=
-void destroy_domain_list(list_t *domain_list)
+We'll be storing the states in a list (see Appendix \ref{sec.lists}),
+so we also define a convenience function for cleaning up.
+<<destroy state list>>=
+void destroy_state_list(list_t *state_list)
{
- domain_list = head(domain_list);
- while (domain_list != NULL)
- destroy_domain((domain_t *) pop(&domain_list));
+ state_list = head(state_list);
+ while (state_list != NULL)
+ destroy_state((state_t *) pop(&state_list));
}
-@
-\subsection{Domain model and group handling}
-
-<<group functions>>=
-<<get tension handler>>
-<<get group>>
-<<int comparison function>>
-<<find possible groups>>
-<<is group active>>
-<<get group list>>
-<<get active groups>>
@
-<<get tension handler>>=
-one_dim_fn_t *get_tension_handler(domain_t *domain)
+We can also define a convenient way to get a state index from it's name.
+<<state index>>=
+int state_index(list_t *state_list, char *name)
{
- if (domain->state == DS_FOLDED)
- return domain->folded_handler;
- else {
- if (domain->state != DS_UNFOLDED)
- fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
- assert(domain->state == DS_UNFOLDED);
- return domain->unfolded_handler;
+ int i=0;
+ state_list = head(state_list);
+ while (state_list != NULL) {
+ if (STRMATCH(S(state_list)->name, name)) return i;
+ state_list = state_list->next;
+ i++;
}
+ return -1;
}
-@
-<<get group>>=
-int get_group(domain_t *domain)
-{
- if (domain->state == DS_FOLDED)
- return domain->folded_group;
- else {
- if (domain->state != DS_UNFOLDED)
- fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
- assert(domain->state == DS_UNFOLDED);
- return domain->unfolded_group;
- }
-}
@
-We already know all possible tension classes, since they are all
-hardcoded into \prog. However, there may be any number of possible
-groups. We define a function [[find_possible_groups]] to search for
-possible (not neccessarily active) groups. It can search for
-subgroups of a single [[handler]], or by setting [[handler]] to
-[[NULL]], subgroups of \emph{any} handler. It returns a list of group
-integers.
-<<find possible groups>>=
-list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
- list_t *ret = NULL;
- list = head(list);
- while (list != NULL) {
- if (handler == NULL || get_tension_handler(D(list)) == handler) {
- push(&ret, &D(list)->folded_group, 1);
- push(&ret, &D(list)->unfolded_group, 1);
- }
- list = list->next;
- }
-
- if (ret == NULL)
- return ret; /* no groups with this handler, no processing needed */
+\subsection{Transition data}
+\label{sec.transition_data}
- /* sort the ret list, and remove duplicates */
- sort(&ret, &int_comparison);
- unique(&ret, &int_comparison);
+Once you've created states (Sections \ref{sec.main},
+\ref{sec.model_selection}, and \ref{sec.state_data}), you need to
+create transitions between them.
+<<transition definitions>>=
+typedef struct transition_struct {
+ state_t *initial_state;
+ state_t *final_state;
+ k_func_t *k; /* transition rate model */
+ void *k_params; /* pointer to k parameters */
+ destroy_data_func_t *destroy_k;
+} transition_t;
+
+/* get the transition data for the current list node */
+#define T(list) ((transition_t *)(list)->d)
+@ [[k]] is a pointer to the function determining the transition rate
+for a given tension. [[k_params]] and [[destroy_k]] have similar
+roles with regards to [[k]] as [[tension_params]] and
+[[destroy_tension]] have with regards to [[tension_handler]] (see
+Section \ref{sec.state_data}).
+
+[[create_]] and [[destroy_transition]] are simple wrappers around
+[[malloc]] and [[free]].
+<<create/destroy transition>>=
+transition_t *create_transition(state_t *initial_state, state_t *final_state,
+ k_func_t *k,
+ void *k_params,
+ destroy_data_func_t *destroy_k)
+{
+ transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
+ assert(ret != NULL);
+ assert(initial_state != NULL);
+ assert(final_state != NULL);
+ ret->initial_state = initial_state;
+ ret->final_state = final_state;
+ ret->k = k;
+ ret->k_params = k_params;
+ ret->destroy_k = destroy_k;
#ifdef DEBUG
- fprintf(stderr, "handler %p possible groups:", handler);
- list = head(ret);
- while (list != NULL) {
- fprintf(stderr, " %d", *((int *)list->d));
- list = list->next;
- }
- fprintf(stderr, "\n");
+ fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
#endif
return ret;
}
-@
-Search a [[list]] of domains, and determine whether a particular model
-class and group number combination is active.
-<<is group active>>=
-int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
+void destroy_transition(transition_t *transition)
{
- list = head(list);
- while (list != NULL) {
- if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
- return 1;
- list = list->next;
+ if (transition) {
+#ifdef DEBUG
+ fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
+#endif
+ if (transition->destroy_k)
+ (*transition->destroy_k)(transition->k_params);
+ free(transition);
}
- return 0;
}
+
@
-Search a [[list]] of domains, and return all domains matching a
-particular model class and group number.
-<<get group list>>=
-list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
+We'll be storing the transitions in a list (see Appendix
+\ref{sec.lists}), so we also define a convenience function for
+cleaning up.
+<<destroy transition list>>=
+void destroy_transition_list(list_t *transition_list)
{
- list_t *ret = NULL;
- list = head(list);
- while (list != NULL) {
- if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
- push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
-#ifdef DEBUG
- fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
-#endif
- }
- list = list->next;
+ transition_list = head(transition_list);
+ while (transition_list != NULL)
+ destroy_transition((transition_t *) pop(&transition_list));
+}
+
+@
+
+\subsection{Graphviz output}
+\label{sec.graphviz_output}
+
+<<print state graph>>=
+void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
+{
+ state_list = head(state_list);
+ transition_list = head(transition_list);
+ fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
+ while (1) {
+ fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
+ if (state_list->next == NULL) break;
+ state_list = state_list->next;
}
- return ret;
+ fprintf(file, "\n");
+ while (transition_list != NULL) {
+ fprintf(file, " n%d -> n%d;\n", state_index(state_list, T(transition_list)->initial_state->name), state_index(state_list, T(transition_list)->final_state->name));
+ transition_list = transition_list->next;
+ }
+ fprintf(file, "}\n");
}
@
-Because all the node data in lists returned by [[get_group_list]] is
-also in the main domain list, you shouldn't destroy the node data
-popped off when destroying the group lists. It will all get cleaned
-up when the main domain list is destroyed.
-Put all this together to scan through a list of domains and construct
-an array of all the active groups.
+\subsection{Domain model and group handling}
+
+<<group functions>>=
+<<crosslink groups>>
+<<get active groups>>
+@
+
+Scan through a list of states and construct an array of all the
+active groups.
<<get active groups>>=
-void get_active_groups(list_t *list,
- int num_tension_handlers, one_dim_fn_t **tension_handlers,
- int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
+void get_active_groups(list_t *state_list,
+ int *pNum_active_groups,
+ one_dim_fn_t ***pPer_group_handlers,
+ void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
{
void **active_groups=NULL;
- one_dim_fn_t *handler, **per_group_handlers=NULL;
- int i, num_active_groups, *group;
- list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
+ one_dim_fn_t *handler, *old_handler, **per_group_handlers=NULL;
+ void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
+ int i, j, num_domains, group, old_group, num_active_groups=0;
+ list_t *active_states_list=NULL;
+ tension_handler_data_t *tdata=NULL;
+ state_t *state;
/* get all the active groups in a list */
- for (i=0; i<num_tension_handlers; i++) {
- handler = tension_handlers[i];
- if (handler == NULL) continue; /* tensionless handler */
- possible_group_numbers = head(find_possible_groups(list, handler));
- while (possible_group_numbers != NULL) {
- group = (int *)pop(&possible_group_numbers);
- if (is_group_active(list, handler, *group) == 1) {
- active_group_list = get_group_list(list, handler, *group);
- push(&active_groups_list, active_group_list, 1);
- push(&per_group_handlers_list, handler, 1);
+ state_list = head(state_list);
#ifdef DEBUG
- fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
- list_print(stderr, active_group_list, "active group");
+ fprintf(stderr, "%s:\t", __FUNCTION__);
+ list_print(stderr, state_list, "states");
#endif
- }
- }
+ while (state_list != NULL) {
+ state = S(state_list);
+ handler = state->tension_handler;
+ group = state->tension_group;
+ num_domains = state->num_domains;
+ if (list_index(active_states_list, state) == -1) {
+ /* we haven't added this state (or it's associated group) yet */
+ if (num_domains > 0 && handler != NULL) { /* add the state */
+ num_active_groups++;
+ push(&active_states_list, state, 1);
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
+#endif
+ for (i=0; i < state->num_grouped_states; i++) {
+ if (state->grouped_states[i]->num_domains > 0) {
+ /* add this related (and active) state now */
+ assert(state->grouped_states[i]->tension_handler == handler);
+ assert(state->grouped_states[i]->tension_group == group);
+ push(&active_states_list, state->grouped_states[i], 1);
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tactive grouped %s state (%p) with %d domain(s)\n", __FUNCTION__, state->grouped_states[i]->name, state->grouped_states[i], state->grouped_states[i]->num_domains);
+#endif
+ }
+ }
+ } /* else empty state or NULL handler */
+ } /* else state already added */
+ state_list = state_list->next;
}
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t", __FUNCTION__);
+ list_print(stderr, active_states_list, "active states");
+#endif
+
+ assert(num_active_groups <= list_length(active_states_list));
+ assert(num_active_groups > 0);
- /* allocate the array we'll be returning */
- num_active_groups = list_length(active_groups_list);
- active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
- assert(active_groups != NULL);
+ /* allocate the arrays we'll be returning */
per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
assert(per_group_handlers != NULL);
+ per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
+ assert(per_group_data != NULL);
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
+#endif
+ for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
+ per_group_data[i] = malloc(sizeof(tension_handler_data_t));
+ assert(per_group_data[i] != NULL);
+#ifdef DEBUG
+ fprintf(stderr, " %d,%p", i, per_group_data[i]);
+#endif
+ }
+#ifdef DEBUG
+ fprintf(stderr, "\n");
+#endif
+
+ /* populate the arrays we'll be returning */
+ active_states_list = head(active_states_list);
+ old_handler = NULL;
+ j = -1; /* j is the current active _group_ index */
+ while (active_states_list != NULL) {
+ state = (state_t *) pop(&active_states_list);
+ handler = state->tension_handler;
+ group = state->tension_group;
+ num_domains = state->num_domains;
+ if (handler != old_handler || group != old_group) {
+ j++; /* move to the next group */
+ tdata = (tension_handler_data_t *) per_group_data[j];
+ per_group_handlers[j] = handler;
+ tdata->group_tension_params = NULL;
+ tdata->env = NULL;
+ tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
+ }
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, group, num_domains);
+#endif
+ for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
+ push(&tdata->group_tension_params, state->tension_params, 1);
+ }
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params);
+#endif
+ old_handler = handler;
+ old_group = group;
+ }
- /* populate the array we'll be returning */
- active_groups_list = head(active_groups_list);
- for (i=0; i<num_active_groups; i++) {
- handler = pop(&per_group_handlers_list);
- per_group_handlers[i] = handler;
-
- active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
- assert(active_groups[i] != NULL);
- active_group_list = pop(&active_groups_list);
- ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
- ((tension_handler_data_t *)active_groups[i])->env = NULL;
- ((tension_handler_data_t *)active_groups[i])->persist = NULL;
+ /* free old groups */
+ if (*pPer_group_handlers != NULL)
+ free(*pPer_group_handlers);
+ if (*pPer_group_data != NULL) {
+ for (i=0; i<*pNum_active_groups; i++)
+ free((*pPer_group_data)[i]);
+ free(*pPer_group_data);
}
- assert(active_groups_list == NULL);
- assert(per_group_handlers_list == NULL);
+ /* set new groups */
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
+#endif
*pNum_active_groups = num_active_groups;
*pPer_group_handlers = per_group_handlers;
- *pActive_groups = active_groups;
+ *pPer_group_data = per_group_data;
}
@
+<<crosslink groups>>=
+void crosslink_groups(list_t *state_groups, list_t *transition_list)
+{
+ list_t *list, *out_trans=NULL, *associates=NULL;
+ one_dim_fn_t *handler;
+ int i, n, group;
+
+ state_groups = head(state_groups);
+ while (state_groups != NULL) {
+ /* find transitions leaving this state */
+#ifdef DEBUG
+ fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
+#endif
+ list = head(transition_list);
+ while (list != NULL) {
+ if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
+ push(&out_trans, T(list), 1);
+ }
+ list = list->next;
+ }
+ n = list_length(out_trans);
+ S(state_groups)->num_out_trans = n;
+ S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
+ assert(S(state_groups)->out_trans != NULL);
+ for (i=0; i<n; i++) {
+ S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
+ }
+
+ /* find states grouped with this state */
+#ifdef DEBUG
+ fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
+#endif
+ handler = S(state_groups)->tension_handler;
+ group = S(state_groups)->tension_group;
+ list = head(state_groups);
+ while (list != NULL) {
+ if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
+ push(&associates, S(list), 1);
+ }
+ list = list->next;
+ }
+ n = list_length(associates);
+ S(state_groups)->num_grouped_states = n;
+ S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
+ assert(S(state_groups)->grouped_states != NULL);
+ for (i=0; i<n; i++) {
+ S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
+ }
+ state_groups = state_groups->next;
+ }
+}
+
+@
\section{String parsing}
-For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
+For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
The model handling in getopt is set up to handle a fixed number of arguments for each model,
so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
need to take care of parsing those parameters themselves.
fail_unless(STRMATCH(param_args[1],"def"), NULL);
}
END_TEST
+START_TEST(test_parse_list_compound)
+{
+ int num_param_args;
+ char **param_args;
+ parse_list_string("abc,{def,ghi}", SEP, '{', '}',
+ &num_param_args, ¶m_args);
+ fail_unless(num_param_args == 2, NULL);
+ fail_unless(STRMATCH(param_args[0],"abc"), NULL);
+ fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
+}
+END_TEST
@
<<parse list string test case defs>>=
TCase *tc_parse_list_string = tcase_create("parse list string");
tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
+tcase_add_test(tc_parse_list_string, test_parse_list_compound);
suite_add_tcase(s, tc_parse_list_string);
@
<<F tests>>
<<determine dt tests>>
<<happens tests>>
-<<does domain unfold tests>>
+<<does domain transition tests>>
<<randomly unfold domains tests>>
<<suite definition>>
@
<<F test case defs>>
<<determine dt test case defs>>
<<happens test case defs>>
- <<does domain unfold test case defs>>
+ <<does domain transition test case defs>>
<<randomly unfold domains test case defs>>
<<F test case add>>
<<determine dt test case add>>
<<happens test case add>>
- <<does domain unfold test case add>>
+ <<does domain transition test case add>>
<<randomly unfold domains test case add>>
/*
Check the searching with [[linear_k]].
Check overwhelming force treatment with the heavyside-esque step [[?]].
+
+Hmm, I'm not really sure what I was doing here...
<<determine dt tests>>=
double linear_k(double F, environment_t *env, void *params)
{
START_TEST(test_determine_dt_linear_k)
{
- environment_t env;
- double dt_max=3.0, Fnot=3.0;
+ state_t folded={0};
+ transition_t unfold={0};
+ environment_t env = {0};
double F[]={0,1,2,3,4,5,6};
- domain_t dom; /* use both parts at once for folded/unfolded */
+ double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
+ list_t *state_list=NULL, *transition_list=NULL;
int i;
env.T = 300.0;
-/*
- dom->next = dom->prev = NULL;
- dom->k_func_t = linear_k;
- dom->folded_params = &Fnot;
- dom->unfolded_params = !!!!!!!!!
- dom->destroy_folded = dom->destroy_unfolded = NULL;
+ folded.tension_handler = &hooke_handler;
+ folded.num_domains = 1;
+ unfold.initial_state = &folded;
+ unfold.k = linear_k;
+ unfold.k_params = &Fnot;
+ push(&state_list, &folded, 1);
+ push(&transition_list, &unfold, 1);
for( i=0; i < sizeof(F)/sizeof(double); i++) {
- fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
+ //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
}
-*/
}
END_TEST
@
<<happens test case add>>=
@
-<<does domain unfold tests>>=
+<<does domain transition tests>>=
@
-<<does domain unfold test case defs>>=
+<<does domain transition test case defs>>=
@
-<<does domain unfold test case add>>=
+<<does domain transition test case add>>=
@
<<randomly unfold domains tests>>=
\section{Balancing group extensions}
-\label{app.tension_balance}
+\label{sec.tension_balance}
For a given total extension $x$ (of the piezo), the various domain
models (WLC, FJC, Hookean springs, \ldots) and groups extend different
<<one dimensional bracket>>
<<one dimensional solve>>
<<x of x0>>
+<<group stiffness function>>
+<<chain stiffness function>>
<<tension balance function>>
@
The entire force balancing problem reduces to a solving two nested
one-dimensional root-finding problems. First we define the one
dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
-non-empty group). $x(x_0)$ is also strictly monotonically increasing,
+populated group). $x(x_0)$ is also strictly monotonically increasing,
so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
stores the last successful value of $x$, since we expect to be taking
small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
<<tension balance function declaration>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
- void **params,
+ void **params, /* array of pointers to tension_handler_data_t */
double last_x,
- double x);
+ double x,
+ double *stiffness);
@
<<tension balance function>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
- void **params,
+ void **params, /* array of pointers to tension_handler_data_t */
double last_x,
- double x)
+ double x,
+ double *stiffness)
{
static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
double F, xo_guess, xo, lb, ub;
assert(params != NULL);
assert(num_tension_groups > 0);
+ if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
+ /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
+ if (x_xo_data.xi != NULL)
+ free(x_xo_data.xi);
+ /* construct new x_xo_data */
+ x_xo_data.n_groups = num_tension_groups;
+ x_xo_data.tension_handler = tension_handlers;
+ x_xo_data.handler_data = params;
+#ifdef DEBUG
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0]);
+#endif
+ x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
+ assert(x_xo_data.xi != NULL);
+ for (i=0; i<num_tension_groups; i++)
+ x_xo_data.xi[i] = last_x;
+ }
+
if (num_tension_groups == 1) { /* only one group, no balancing required */
xo = x;
+ x_xo_data.xi[0] = xo;
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
+#endif
} else {
- //fprintf(stderr, "balancing for x = %g with ", x);
- //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
- //fprintf(stderr, "\n");
- if (last_x == -1) { /* new group setup, reset x_xo_data */
- /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
- if (x_xo_data.xi != NULL)
- free(x_xo_data.xi);
- /* construct new x_xo_data */
- x_xo_data.n_groups = num_tension_groups;
- x_xo_data.tension_handler = tension_handlers;
- x_xo_data.handler_data = params;
- x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
- for (i=0; i<num_tension_groups; i++)
- x_xo_data.xi[i] = -1.0;
- }
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
+ for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
+ fprintf(stderr, "\n");
+#endif
if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
if (x == 0) xo_guess = length_scale;
else xo_guess = x/num_tension_groups;
}
#endif
}
+
F = (*tension_handlers[0])(xo, params[0]);
+
+ if (stiffness != NULL)
+ *stiffness = chain_stiffness(&x_xo_data, min_relx);
+
return F;
}
+
@
<<tension balance internal definitions>>=
<<x of x0 definitions>>=
typedef struct x_of_xo_data_struct {
int n_groups;
- one_dim_fn_t **tension_handler; /* array of fn pointers */
- void **handler_data; /* array of void* pointers */
- double *xi;
+ one_dim_fn_t **tension_handler; /* array of fn pointers */
+ void **handler_data; /* array of pointers to tension_handler_data_t */
+ double *xi; /* array of group extensions */
} x_of_xo_data_t;
@
double x_of_xo(double xo, void *pdata)
{
x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
- double F, x=0, xi, xi_guess, lb, ub;
+ double F, x=0, xi, xi_guess, lb, ub, slope;
int i;
double min_relx=1e-6, min_rely=1e-6;
int max_steps=200;
assert(data->n_groups > 0);
data->xi[0] = xo;
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0]);
+ fflush(stderr);
+#endif
F = (*data->tension_handler[0])(xo, data->handler_data[0]);
#ifdef DEBUG
- fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
+ fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
+ fflush(stderr);
#endif
x += xo;
if (data->xi) data->xi[0] = xo;
if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
else xi_guess = data->xi[i];
#ifdef DEBUG
- fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
+ fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
#endif
oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
#ifdef DEBUG
- fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
+ fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
#endif
data->xi[i] = xi;
x += xi;
if (data->xi) data->xi[i] = xi;
}
#ifdef DEBUG
- fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
+ fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
#endif
return x;
}
@
+<<group stiffness function>>=
+double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
+{
+ double dx, xl, F, Fl, stiffness;
+ assert(x >= 0);
+ assert(relx > 0 && relx < 1);
+ if (x == 0 || relx == 0) {
+ dx = 1e-12; /* HACK, how to get length scale? */
+ xl = x-dx;
+ if (xl < 0) {
+ xl = 0;
+ x = xl+dx;
+ }
+ } else {
+ dx = x*relx;
+ xl = x-dx;
+ }
+ F = (*tension_handler)(x, handler_data);
+ Fl = (*tension_handler)(xl, handler_data);
+ stiffness = (F-Fl)/dx;
+ assert(stiffness > 0);
+ return stiffness;
+}
+@
+
+<<chain stiffness function>>=
+double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
+{
+ double stiffness;
+ int i;
+ /* add group stiffnesses in series */
+ for (i=0; i < x_xo_data->n_groups; i++) {
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
+ fflush(stderr);
+#endif
+ stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
+ }
+ stiffness = 1.0 / stiffness;
+ return stiffness;
+}
+
+@
+
Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
number of steps for monotonically increasing $f(x)$. Simple
bisection, so it's robust and converges fairly quickly. You can set
/* check some simple cases */
if (ylb == yub) {
- if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
+ if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
else return lb; /* any x on the range [lb, ub] would work */
}
if (ylb == y) { x=lb; yx=ylb; goto end; }
double k=5, x=3, last_x=2.0, F;
one_dim_fn_t *handlers[] = {&hooke};
void *data[] = {&k};
- F = tension_balance(1, handlers, data, last_x, x);
+ F = tension_balance(1, handlers, data, last_x, x, NULL);
fail_unless(F = k*x, NULL);
}
END_TEST
@
We can also test balancing two springs with different spring constants.
-The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
+The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
<<tension balance function tests>>=
START_TEST(test_double_hooke)
{
double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
one_dim_fn_t *handlers[] = {&hooke, &hooke};
void *data[] = {&k1, &k2};
- F = tension_balance(2, handlers, data, last_x, x);
+ F = tension_balance(2, handlers, data, last_x, x, NULL);
x1e = x*k2/(k1+k2);
Fe = k1*x1e;
x2e = x1e*k1/k2;
@
\section{Lists}
+\label{sec.lists}
The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
<<list to array declaration>>
<<move declaration>>
<<sort declaration>>
+<<list index declaration>>
+<<list element declaration>>
<<unique declaration>>
<<list print declaration>>
@
<<list to array>>
<<move>>
<<sort>>
+<<list index>>
+<<list element>>
<<unique>>
<<list print>>
@
else if (a == b) return 0;
else return -1;
}
+
@
<<sort declaration>>=
void sort(list_t **pList, comparison_fn_t *comp);
}
*pList = list;
}
+
@
+[[list_index]] finds the location of [[data]] in [[list]]. Returns
+[[-1]] if [[data]] is not in [[list]].
+<<list index declaration>>=
+int list_index(list_t *list, void *data);
+@
+
+<<list index>>=
+int list_index(list_t *list, void *data)
+{
+ int i=0;
+ list = head(list);
+ while (list != NULL) {
+ if (list->d == data) return i;
+ list = list->next;
+ i++;
+ }
+ return -1;
+}
+
+@
+
+[[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
+<<list element declaration>>=
+void *list_element(list_t *list, int i);
+@
+<<list element>>=
+void *list_element(list_t *list, int i)
+{
+ int j=0;
+ list = head(list);
+ while (list != NULL) {
+ if (i == j) return list->d;
+ list = list->next;
+ j++;
+ }
+ assert(0==1);
+}
+
+@
+
+
[[unique]] removes duplicates from a sorted list according to a user
specified comparison. Don't do this unless you have other pointers
any dynamically allocated data in your list, because [[unique]] just
Most of this is also POSIX-specific (as far as I know), so you'll have to do some legwork to port it to non-POSIX systems.
We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
-If you feel the benefits do \emph{not} outweigh the costs, or you are on a non-POSIX system (e.g. MS Windows without Cygwin), you should probably hardcode your functions in \lang.
+If you feel the benefits do \emph{not} outweigh the costs, or you are on a non-POSIX system (e.g.\ MS Windows without Cygwin), you should probably hardcode your functions in \lang.
Then you can either statically or dynamically link to those hardcoded functions.
While much less flexible, this approach would be a fairly simple-to-implement fallback.
The [[alarm]] calls set and clear a timeout on the returned output.
If the timeout expires, the process would get a [[SIGALRM]], but it doesn't have a [[SIGALRM]] handler, so it gets a [[SIGKILL]] and dies.
This protects against invalid input for which a line of output is not printed to [[stdout]].
-Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
+Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
If you are getting strange results, check your python code seperately. TODO, better error handling.
Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
double const_tension_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double F;
assert (x >= 0.0);
$$
F = \p({\sum_i \frac{1}{k_i}})^{-1} x
$$
-For a simple proof, see Appendix \ref{app.math_hooke}.
+For a simple proof, see Appendix \ref{sec.math_hooke}.
<<hooke tension function declaration>>=
extern double hooke_handler(double x, void *pdata);
double hooke_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double k=0.0;
assert (x >= 0.0);
}
k = 1.0 / k;
#ifdef DEBUG
- fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
- __FUNCTION__, k, x, k*x);
+ fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
+ __FUNCTION__, k, x, k*x, data);
#endif
return k*x;
}
double wlc_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double p, L=0.0;
list = head(list);
double fjc_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double l;
int N=0;
printf("#initializing model %s with parameters %s\n", model->name, model->params);
}
INIT_MODEL("utility", model, model->params, params);
- tdata.group = NULL;
- push(&tdata.group, params, 1);
+ tdata.group_tension_params = NULL;
+ push(&tdata.group_tension_params, params, 1);
tdata.env = &env;
tdata.persist = NULL;
if (model->handler == NULL) {
printf("Reached large force limit %g > %g\n", F, Fmax);
}
}
- params = pop(&tdata.group);
+ params = pop(&tdata.group_tension_params);
if (model->destructor)
(*model->destructor)(params);
return 0;
bell_param_t *p = (bell_param_t *) bell_params;
assert(F >= 0); assert(env->T > 0);
assert(p != NULL);
- assert(p->knot > 0); assert(p->dx >= 0);
+ assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
/*
printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
doesn't hold for biotin-streptavadin binding, and I extended their
results to I27 for the 2009 Biophysical Society Annual
Meeting\cite{walton08,king09}. More details to follow when I get done
-with the conference...
+with the conference\ldots
We adjust Bell's model to
$$
Which we initialized with parameters for ddFLN4\citep{schlierf06}.
There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
-Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{app.adaptive_dt} reduces it's timestep appropriately.
+Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{sec.adaptive_dt} reduces it's timestep appropriately.
The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
It works so long as [[val_1]] is not [[false]].
\section{Math}
\subsection{Hookean springs in series}
-\label{app.math_hooke}
+\label{sec.math_hooke}
\begin{theorem}
The effective spring constant for $n$ springs $k_i$ in series is given by
\bibliography{sawsim}
\end{document}
-
projects / sawsim.git / commitdiff