%%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
-% we have our own preamble,
+% we have our own preamble,
% use `noweave -delay` to not print noweb's automatic one
% -*- mode: Noweb; noweb-code-mode: c-mode -*-
\documentclass[letterpaper, 11pt]{article}
tables to accelerate slow unfolding rate model evaluation, and command
line control of tension groups.
+Version 0.5 added the stiffness environmental parameter and it's
+associated unfolding models.
+
<<version definition>>=
-#define VERSION "0.4"
-@
+#define VERSION "0.5"
+@
\subsection{License}
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
The author may be contacted at <wking@drexel.edu> on the Internet, or
write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
Philadelphia PA 19104, USA.
-@
+@
<<license comment>>=
/*
<<license>>
*/
-@
+@
\section{Main}
\begin{packed_item}
\item Null (Appendix \ref{sec.null_k}),
\item Bell's model (Appendix \ref{sec.bell}),
+ \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
\item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
\item Kramers' saddle point model (Appendix \ref{sec.kramers}).
\end{packed_item}
get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
setup();
- if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
+ 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);
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);
- if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
+ if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
if (xstep == 0) {
x += v*dt;
} else {
Environmental parameters important in determining reaction rates and
tensions (e.g. temperature, pH) are stored in a single structure to
facilitate extension to more complicated models in the future.
-<<environment definition>>=
+<<environment definition>>=
typedef struct environment_struct {
double T;
+ double stiffness;
} environment_t;
-@
+@
<<global.h>>=
#ifndef GLOBAL_H
<<create/destroy definitions>>
<<model globals>>
#endif /* GLOBAL_H */
-@
+@
Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
included multiple times.
<<happens>>
<<does domain unfold>>
<<randomly unfold domains>>
-@
+@
\subsection{Tension}
\label{sec.find_tension}
<<find tension definitions>>=
#define NUM_TENSION_MODELS 5
-@
+@
<<tension handler array global declaration>>=
extern one_dim_fn_t *tension_handlers[];
@
-<<tension handler array global>>=
+<<tension handler array global>>=
one_dim_fn_t *tension_handlers[] = {
NULL,
&const_tension_handler,
&fjc_handler,
};
-@
+@
<<tension model global declarations>>=
<<tension handler array global declaration>>
-@
+@
<<tension handler types>>=
typedef struct tension_handler_data_struct {
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
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]].
+[[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
+\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)
{
list_print(stderr, domain_list, "domain list");
#endif
- if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
+ if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
/* free old statics */
if (per_group_handlers != NULL)
free(per_group_handlers);
} /* else roll 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
+ }
+
last_x = x;
+
return F;
}
@ For the moment, we will restrict our group boundaries to a single
shouldn't be a problem. To reassure yourself, you can ask the
simulator to print the smallest timestep that was used with TODO.
<<randomly unfold domains>>=
-int random_unfoldings(list_t *domain_list, double tension,
+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
+ 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);
}
return 0;
}
-@
+@
\subsection{Adaptive timesteps}
\label{sec.adaptive_dt}
\ref{app.adaptive_dt}.
\section{Models}
-
+
TODO: model intro.
The models provide the physics of the simulation, but the simulation
<<model globals>>=
#define k_B 1.3806503e-23 /* J/K */
-@
+@
\section{Command line interface}
<<index k model>>
<<generate domain>>
<<get options>>
-@
+@
\subsection{Model selection}
\label{app.model_selection}
<<get options definitions>>=
<<model getopt definitions>>
-@
+@
<<create/destroy definitions>>=
typedef void *create_data_func_t(char **param_strings);
typedef void destroy_data_func_t(void *param_struct);
-@
+@
<<model getopt definitions>>=
<<tension model getopt definitions>>
<<k model getopt definitions>>
-@
+@
\subsubsection{Tension}
create_data_func_t *creator; /* param-string -> model-data-struct fn */
destroy_data_func_t *destructor; /* fn to free the model data structure */
} tension_model_getopt_t;
-@
+@
<<tension model getopt array definition>>=
tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
<<worm-like chain tension model getopt>>,
<<freely-jointed chain tension model getopt>>
};
-@
+@
\subsubsection{Unfolding rate}
<<k model getopt definitions>>=
-#define NUM_K_MODELS 5
+#define NUM_K_MODELS 6
typedef struct k_model_getopt_struct {
char *name;
create_data_func_t *creator;
destroy_data_func_t *destructor;
} k_model_getopt_t;
-@
+@
<<k model getopt array definition>>=
k_model_getopt_t k_models[NUM_K_MODELS] = {
<<null k model getopt>>,
<<const k model getopt>>,
<<bell k model getopt>>,
+ <<kbell k model getopt>>,
<<kramers k model getopt>>,
<<kramers integ k model getopt>>
};
-@
+@
\subsection{help}
printf(" 123.456e-12\n");
printf(" ...\n");
printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
- printf(" #Position (m)\tForce (N)\n");
+ printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
printf(" 0.001\t0.0023\n");
printf(" ...\n");
printf("-V\tChange output to verbose mode\n");
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);
}
-@
+@
\subsection{Get options}
assert(env->T > 0.0);
return;
}
-@
+@
<<index tension model>>=
int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
}
return i;
}
-@
+@
<<index k model>>=
int index_k_model(int n_models, k_model_getopt_t *models, char *name)
}
return i;
}
-@
+@
<<initialize model definition>>=
/* requires int num_param_args and char **param_args in the current scope
param_pointer = NULL; \
free_parsed_list(num_param_args, param_args); \
} while (0);
-@
+@
<<generate domain>>=
void *generate_domain(enum domain_state_t state,
unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
);
}
-@
+@
\phantomsection
\appendix
<<globals>>
<<functions>>
<<main program>>
-@
+@
We include [[math.h]], so don't forget to link to the libm with `-lm'.
<<includes>>=
#include "tension_model.h"
#include "parse.h"
#include "accel_k.h"
-@
+@
<<definitions>>=
<<version definition>>
<<initialize model definition>>
<<get options definitions>>
<<domain definitions>>
-@
+@
<<globals>>=
<<flag globals>>
<<model globals>>
-@
+@
<<functions>>=
<<create/destroy domain>>
<<group functions>>
<<simulation functions>>
<<boilerplate functions>>
-@
+@
<<boilerplate functions>>=
<<setup>>
<<get option functions>>
-@
+@
<<setup>>=
void setup(void)
{
srand(getpid()*time(NULL)); /* seed rand() */
}
-@
+@
<<flag definitions>>=
/* in octal b/c of prefixed '0' */
<<max/min definitions>>=
#define MAX(a,b) ((a)>(b) ? (a) : (b))
#define MIN(a,b) ((a)<(b) ? (a) : (b))
-@
+@
Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
<<string comparison definition and globals>>=
strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
!strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
/* +1 to also compare the '\0' */
-@
+@
We also define a macro for our [[check]] unit testing
<<check relative error>>=
-(received-expected)/expected, max_err, #received, \
expected, received); \
} while(0)
-@
+@
<<happens>>=
int happens(double probability)
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,
+ list_t *domain_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.
//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");
+ //printf("overshot max_dt\n");
dtU = max_dt;
}
if (v == 0) /* discrete stepping, so force is temporatily constant */
//printf("Try for dt = %g (F = %g)\n", dt, F);
k = accel_k(D(domain_list)->k, F, env, D(domain_list)->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);
+ //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);
//printf("Try for dt = %g (F = %g)\n", dt, F);
k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
- //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
+ //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 */
dtUCur = target_prob / k; /* safe timestep back from x+dtU */
if (dtCur > dt) /* safe timestep back from x+dt covers dt */
dtL = dt;
else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
- dtU = dt; /* ... stepping out only dtCur would be safe */
+ dtU = dt; /* ... stepping out only dtCur would be safe */
dtUCur = dtCur;
} else
break; /* dtCur = 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.
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);
}
return dt;
}
-@
+@
\subsection{Domain data}
#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.
<<is group active>>
<<get group list>>
<<get active groups>>
-@
+@
<<get tension handler>>=
one_dim_fn_t *get_tension_handler(domain_t *domain)
return domain->unfolded_handler;
}
}
-@
+@
<<get group>>=
int get_group(domain_t *domain)
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
#endif
return ret;
}
-@
+@
Search a [[list]] of domains, and determine whether a particular model
class and group number combination is active.
}
return 0;
}
-@
+@
Search a [[list]] of domains, and return all domains matching a
particular model class and group number.
}
return ret;
}
-@
+@
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
}
}
- /* allocate the array we'll be returning */
+ /* 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);
*pPer_group_handlers = per_group_handlers;
*pActive_groups = active_groups;
}
-@
+@
\section{String parsing}
<<parse definitions>>
<<parse declarations>>
#endif /* PARSE */
-@
+@
<<parse module makefile lines>>=
NW_SAWSIM_MODS += parse
CHECK_BINS += check_parse
-check_parse_MODS =
-@
+check_parse_MODS =
+@
<<parse definitions>>=
#define SEP ',' /* argument separator character */
-@
+@
<<parse declarations>>=
extern void parse_list_string(char *string, char sep, char deeper, char shallower,
int *num, char ***string_array);
extern void free_parsed_list(int num, char **string_array);
-@
+@
[[parse_list_string]] allocates memory, don't forget to free it
afterward with [[free_parsed_list]]. It does not alter the original.
<<parse includes>>
#include "parse.h"
<<parse delimited list string functions>>
-@
+@
<<parse includes>>=
#include <assert.h> /* assert() */
#include <stdio.h> /* fprintf(), stdout *//*!!*/
#include <string.h> /* strlen() */
#include "parse.h"
-@
+@
\subsection{Parsing unit tests}
<<check relative error>>
<<parse test suite>>
<<main check program>>
-@
+@
<<parse check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
#include <string.h> /* strlen() */
<<check includes>>
#include "parse.h"
-@
+@
<<parse test suite>>=
<<parse list string tests>>
<<parse suite definition>>
-@
+@
<<parse suite definition>>=
Suite *test_suite (void)
<<parse list string test case add>>
return s;
}
-@
+@
<<parse list string tests>>=
<<functions>>
<<test suite>>
<<main check program>>
-@
+@
<<check includes>>=
#include <check.h>
-@
+@
<<check globals>>=
-@
+@
<<test suite>>=
<<F tests>>
<<does domain unfold tests>>
<<randomly unfold domains tests>>
<<suite definition>>
-@
+@
<<suite definition>>=
Suite *test_suite (void)
*/
return s;
}
-@
+@
<<main check program>>=
int main(void)
srunner_free(sr);
return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
-@
+@
\subsection{F tests}
<<F tests>>=
<<worm-like chain tests>>
-@
+@
<<F test case defs>>=
<<worm-like chain test case def>>
-@
+@
<<F test case add>>=
<<worm-like chain test case add>>
-@
+@
<<worm-like chain tests>>=
extern double wlc(double x, double T, double p, double L);
* = 0.25 * 3 = 3/4
* linear term = x/L = 0.5
* nonlinear + linear = 0.75 + 0.5 = 1.25
- * wlc = 10e21*1.25 = 12.5e21
+ * wlc = 10e21*1.25 = 12.5e21
*/
fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
"wlc(%g, %g, %g, %g) = %g != %g",
x, T, p, L, wlc(x, T, p, L), 12.5e21);
}
END_TEST
-@
+@
<<worm-like chain test case def>>=
TCase *tc_wlc = tcase_create("WLC");
-@
+@
<<worm-like chain test case add>>=
tcase_add_test(tc_wlc, test_wlc_at_zero);
tcase_add_test(tc_wlc, test_wlc_at_half);
suite_add_tcase(s, tc_wlc);
-@
+@
\subsection{Model tests}
<<license comment>>
<<tension balance function declaration>>
#endif /* TENSION_BALANCE_H */
-@
+@
<<tension balance functions>>=
<<one dimensional bracket>>
<<one dimensional solve>>
<<x of x0>>
<<tension balance function>>
-@
+@
<<tension balance module makefile lines>>=
NW_SAWSIM_MODS += tension_balance
CHECK_BINS += check_tension_balance
-check_tension_balance_MODS =
-@
+check_tension_balance_MODS =
+@
The entire force balancing problem reduces to a solving two nested
one-dimensional root-finding problems. First we define the one
void **params,
double last_x,
double x);
-@
+@
<<tension balance function>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
F = (*tension_handlers[0])(xo, params[0]);
return F;
}
-@
+@
<<tension balance internal definitions>>=
<<x of x0 definitions>>
-@
+@
<<x of x0 definitions>>=
typedef struct x_of_xo_data_struct {
void **handler_data; /* array of void* pointers */
double *xi;
} x_of_xo_data_t;
-@
+@
<<x of x0>>=
double x_of_xo(double xo, void *pdata)
#endif
return x;
}
-@
+@
Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
number of steps for monotonically increasing $f(x)$. Simple
<<one dimensional function definition>>=
/* equivalent to gsl_func_t */
typedef double one_dim_fn_t(double x, void *params);
-@
+@
<<one dimensional solve declaration>>=
double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
- double min_relx, double min_rely, int max_steps,
+ double min_relx, double min_rely, int max_steps,
double *pYx);
-@
+@
<<one dimensional solve>>=
/* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
- double min_relx, double min_rely, int max_steps,
+ double min_relx, double min_rely, int max_steps,
double *pYx)
{
double yx, ylb, yub, x;
assert(ub >= lb);
ylb = (*f)(lb, params);
yub = (*f)(ub, params);
-
+
/* check some simple cases */
if (ylb == yub) {
if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
#endif
return x;
}
-@
+@
The one dimensional solver needs a bracketed solution, which sometimes we don't have.
Generate bracketing $x$ values through bisection or exponential growth.
<<one dimensional bracket declaration>>=
void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
-@
+@
<<one dimensional bracket>>=
void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
}
//fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
}
-@
+@
\subsection{Balancing implementation}
<<tension balance internal definitions>>
<<tension balance functions>>
-@
+@
<<tension balance includes>>=
#include <assert.h> /* assert() */
#include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
#include <stdio.h> /* fprintf(), stdout */
#include "global.h"
-@
+@
\subsection{Balancing unit tests}
<<tension balance check includes>>
<<tension balance test suite>>
<<main check program>>
-@
+@
<<tension balance check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
<<check includes>>
#include "global.h"
#include "tension_balance.h"
-@
+@
<<tension balance test suite>>=
<<tension balance function tests>>
<<tension balance suite definition>>
-@
+@
<<tension balance suite definition>>=
Suite *test_suite (void)
<<tension balance function test case adds>>
return s;
}
-@
+@
<<tension balance function tests>>=
<<check relative error>>
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}.
CHECK_ERR(1e-6, Fe, F);
}
END_TEST
-@
+@
TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
<<tension balance function test case defs>>=
TCase *tc_tbfunc = tcase_create("tension balance function");
-@
+@
<<tension balance function test case adds>>=
tcase_add_test(tc_tbfunc, test_single_function);
tcase_add_test(tc_tbfunc, test_double_hooke);
suite_add_tcase(s, tc_tbfunc);
-@
+@
\section{Lists}
<<list definitions>>
<<list declarations>>
#endif /* LIST_H */
-@
+@
<<list declarations>>=
<<head and tail declarations>>
<<sort declaration>>
<<unique declaration>>
<<list print declaration>>
-@
+@
<<list functions>>=
<<create and destroy node>>
<<sort>>
<<unique>>
<<list print>>
-@
+@
<<list module makefile lines>>=
NW_SAWSIM_MODS += list
CHECK_BINS += check_list
-check_list_MODS =
-@
+check_list_MODS =
+@
<<list definitions>>=
typedef struct list_struct {
} list_t;
<<comparison function typedef>>
-@
+@
[[head]] and [[tail]] return pointers to the head and tail nodes of the list:
<<head and tail declarations>>=
list_t *head(list_t *list);
list_t *tail(list_t *list);
-@
+@
<<head and tail>>=
list_t *head(list_t *list)
{
}
return list;
}
-@
+@
<<list length declaration>>=
int list_length(list_t *list);
-@
+@
<<list length>>=
int list_length(list_t *list)
{
}
return i;
}
-@
+@
[[push]] inserts a node relative to the current position in the list
without changing the current position.
otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
<<push declaration>>=
void push(list_t **pList, void *data, int below);
-@
+@
<<push>>=
void push(list_t **pList, void *data, int below)
-{
+{
list_t *list, *new_node;
assert(pList != NULL);
list = *pList;
new_node->prev = list;
}
}
-@
+@
[[pop]] removes the current domain node, moving the current position
to the node after it, unless that node is [[NULL]], in which case move
the current position to the node before it.
<<pop declaration>>=
void *pop(list_t **pList);
-@
+@
<<pop>>=
void *pop(list_t **pList)
{
destroy_node(popped);
return data;
}
-@
+@
[[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
If the cleanup function is [[NULL]], no cleanup function is called.
The nodes are not popped in any particular order.
<<list destroy declaration>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy);
-@
+@
<<list destroy>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy)
{
destroy(data);
}
}
-@
+@
[[list_to_array]] converts a list to an array of pointers, preserving order.
<<list to array declaration>>=
void list_to_array(list_t *list, int *length, void ***array);
-@
+@
<<list to array>>=
void list_to_array(list_t *list, int *pLength, void ***pArray)
{
*pLength = length;
*pArray = array;
}
-@
+@
[[move]] moves the current element along the list in either direction.
<<move declaration>>=
void move(list_t **pList, int downstream);
-@
+@
<<move>>=
void move(list_t **pList, int downstream)
{
list = flipper;
if (downstream == 0) {
push(&list, data, 0); /* b>A>c>d */
- list = list->prev; /* B>a>c>d */
+ list = list->prev; /* B>a>c>d */
} else {
push(&list, data, 1); /* a>C>b>d */
list = list->next; /* a>c>B>d */
}
*pList = list;
}
-@
+@
[[sort]] sorts a list from smallest to largest according to a user
specified comparison.
<<comparison function typedef>>=
typedef int (comparison_fn_t)(void *A, void *B);
-@
+@
For example
<<int comparison function>>=
int int_comparison(void *A, void *B)
else if (a == b) return 0;
else return -1;
}
-@
+@
<<sort declaration>>=
void sort(list_t **pList, comparison_fn_t *comp);
-@
+@
Here we use the bubble sort algorithm.
<<sort>>=
void sort(list_t **pList, comparison_fn_t *comp)
}
*pList = list;
}
-@
+@
[[unique]] removes duplicates from a sorted list according to a user
specified comparison. Don't do this unless you have other pointers
drops any non-unique data without freeing it.
<<unique declaration>>=
void unique(list_t **pList, comparison_fn_t *comp);
-@
+@
<<unique>>=
void unique(list_t **pList, comparison_fn_t *comp)
{
}
*pList = list;
}
-@
+@
[[list_print]] prints a list to a [[FILE *]] stream.
<<list print declaration>>=
void list_print(FILE *file, list_t *list, char *list_name);
-@
+@
<<list print>>=
void list_print(FILE *file, list_t *list, char *list_name)
{
}
fprintf(file, "\n");
}
-@
+@
[[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
<<create and destroy node>>=
<<list includes>>
#include "list.h"
<<list functions>>
-@
+@
<<list includes>>=
#include <assert.h> /* assert() */
#include <stdlib.h> /* malloc(), free(), rand() */
#include <stdio.h> /* fprintf(), stdout, FILE */
#include "global.h" /* destroy_data_func_t */
-@
+@
\subsection{List unit tests}
<<list check includes>>
<<list test suite>>
<<main check program>>
-@
+@
<<list check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
<<check includes>>
#include "global.h"
#include "list.h"
-@
+@
<<list test suite>>=
<<push tests>>
<<pop tests>>
<<list suite definition>>
-@
+@
<<list suite definition>>=
Suite *test_suite (void)
<<pop test case adds>>
return s;
}
-@
+@
<<push tests>>=
START_TEST(test_push)
}
}
END_TEST
-@
+@
<<push test case defs>>=
TCase *tc_push = tcase_create("push");
-@
+@
<<push test case adds>>=
tcase_add_test(tc_push, test_push);
suite_add_tcase(s, tc_push);
-@
+@
<<pop tests>>=
-@
+@
<<pop test case defs>>=
-@
+@
<<pop test case adds>>=
-@
+@
\section{Function string evaluation}
<<string eval function declaration>>
<<string eval teardown declaration>>
#endif /* STRING_EVAL_H */
-@
+@
<<string eval module makefile lines>>=
NW_SAWSIM_MODS += string_eval
CHECK_BINS += check_string_eval
-check_string_eval_MODS =
-@
+check_string_eval_MODS =
+@
For an introduction to POSIX process control, see\\
\url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
//#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
//static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
-@
+@
The [[i]] option lets Python know that it should run in interactive mode.
In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
In interactive mode, python acts on every instruction as soon as it is recieved.
-The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it.
+The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it.
%The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply [[pass]]es so that the silly Python header information is not printed. We leave the prompts in, because we scan for them to determine when the output has completed.
Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
#define STDIN 0 /* initial index of stdin pair */
#define STDOUT 2 /* initial index of stdout pair */
-@
+@
So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
As a finishing touch, we can promote the POSIX file descriptors ([[read]]/[[write]] interface) into the more familiar [[stdio.h]] \emph{streams} ([[fprintf]]/[[fgetc]] interface) using [[fdopen]], which creates a stream from an open file descriptor.
<<string eval setup declaration>>=
extern void string_eval_setup(FILE **pIn, FILE **pOut);
-@
+@
<<string eval setup definition>>=
void string_eval_setup(FILE **pIn, FILE **pOut)
{
int rc;
assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
-
+
pid = fork(); /* split process into two copies */
if (pid == -1) { /* fork error */
As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
<<string eval teardown declaration>>=
extern void string_eval_teardown(FILE **pIn, FILE **pOut);
-@
+@
<<string eval teardown definition>>=
void string_eval_teardown(FILE **pIn, FILE **pOut)
}
*/
}
-@
+@
The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
\subsection{String evaluation implementation}
#include "string_eval.h"
<<string eval internal definitions>>
<<string eval functions>>
-@
+@
<<string eval includes>>=
#include <assert.h> /* assert() */
#include <sys/types.h> /* pid_t */
#include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
#include <wait.h> /* wait() */
-@
+@
<<string eval internal definitions>>=
<<string eval subprocess definitions>>
<<string eval pipe definitions>>
-@
+@
<<string eval functions>>=
<<string eval setup definition>>
<<string eval function definition>>
<<string eval teardown definition>>
-@
+@
\subsection{String evaluation unit tests}
<<string comparison definition and globals>>
<<string eval test suite>>
<<main check program>>
-@
+@
<<string eval check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
#include <signal.h> /* SIGKILL */
<<check includes>>
#include "string_eval.h"
-@
+@
<<string eval test suite>>=
<<string eval tests>>
<<string eval suite definition>>
-@
+@
<<string eval suite definition>>=
Suite *test_suite (void)
<<string eval test case add>>
return s;
}
-@
+@
<<string eval tests>>=
START_TEST(test_setup_teardown)
char input[80], output[80]={};
string_eval_setup(&in, &out);
sprintf(input, "print ABCDefg\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
string_eval_teardown(&in, &out);
}
END_TEST
char input[80], output[80]={};
string_eval_setup(&in, &out);
sprintf(input, "print 3+4*5\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"23\n"), NULL);
string_eval_teardown(&in, &out);
}
char input[80], output[80]={};
string_eval_setup(&in, &out);
sprintf(input, "print 3+4*5\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"23\n"), NULL);
sprintf(input, "print (3**2 + 4**2)**0.5\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"5.0\n"), NULL);
string_eval_teardown(&in, &out);
}
sprintf(input, "print 3+4*5\n");
fprintf(in, "A = 3\n");
sprintf(input, "print A*3\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"9\n"), NULL);
string_eval_teardown(&in, &out);
}
\section{Accelerating function evaluation}
My first version-0.3 code was running very slowly.
-With the options suggested in the help
-([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
+With the options suggested in the help
+([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
That's only 15 calls per solution, so the search algorithm seems reasonable.
double accel_k(k_func_t *k, double F, environment_t *env, void *params);
void free_accels();
#endif /* ACCEL_K_H */
-@
+@
<<accel k module makefile lines>>=
NW_SAWSIM_MODS += accel_k
#CHECK_BINS += check_accel_k
-check_accel_k_MODS =
-@
+check_accel_k_MODS =
+@
<<accel-k.c>>=
#include <assert.h> /* assert() */
{
int i=num_accels;
accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
- assert(accels != NULL);
+ assert(accels != NULL);
accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
accels[i].k = k;
accels[i].env = env;
* WARNING: we're assuming param and env are the same. */
return interp_table_eval(accels[i].itable, accels+i, F);
}
- }
+ }
/* set up a new acceleration environment */
i = add_accel_k(k, env, params);
return interp_table_eval(accels[i].itable, accels+i, F);
}
-@
+@
\section{Tension models}
\label{sec.tension_models}
<<find tension definitions>>
<<tension model global declarations>>
#endif /* TENSION_MODEL_H */
-@
+@
<<tension model module makefile lines>>=
NW_SAWSIM_MODS += tension_model
#CHECK_BINS += check_tension_model
-check_tension_model_MODS =
-@
+check_tension_model_MODS =
+@
<<tension model utils makefile lines>>=
TENSION_MODEL_MODS = tension_model parse list tension_balance
TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
<<null tension model getopt>>=
{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
-@
+@
\subsection{Constant}
\label{sec.const_tension}
<<constant tension functions>>=
<<constant tension function>>
<<constant tension parameter create/destroy functions>>
-@
+@
<<constant tension model declarations>>=
<<constant tension function declaration>>
<<constant tension parameter create/destroy function declarations>>
<<constant tension model global declarations>>
-@
+@
An infinitely stretchable domain providing a constant tension.
<<constant tension function declaration>>=
extern double const_tension_handler(double x, void *pdata);
-@
+@
<<constant tension function>>=
double const_tension_handler(double x, void *pdata)
{
return F;
}
-@
+@
<<constant tension structure>>=
typedef struct const_tension_param_struct {
double F; /* tension (force) in N */
} const_tension_param_t;
-@
+@
<<constant tension parameter create/destroy function declarations>>=
extern void *string_create_const_tension_param_t(char **param_strings);
extern void destroy_const_tension_param_t(void *p);
-@
+@
<<constant tension parameter create/destroy functions>>=
const_tension_param_t *create_const_tension_param_t(double F)
extern int num_const_tension_params;
extern char *const_tension_param_descriptions[];
extern char const_tension_param_string[];
-@
+@
<<constant tension model globals>>=
int num_const_tension_params = 1;
char *const_tension_param_descriptions[] = {"tension F, N"};
char const_tension_param_string[] = "0";
-@
+@
<<constant tension model getopt>>=
{&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
-@
+@
\subsection{Hooke}
\label{sec.hooke}
<<hooke functions>>=
<<hooke function>>
<<hooke parameter create/destroy functions>>
-@
+@
<<hooke tension model declarations>>=
<<hooke tension function declaration>>
<<hooke tension parameter create/destroy function declarations>>
<<hooke tension model global declarations>>
-@
+@
The tension of a single spring is given by $F=kx$ for some spring constant $k$.
The behavior of a series of springs $k_i$ in series is given by
<<hooke tension function declaration>>=
extern double hooke_handler(double x, void *pdata);
-@
+@
<<hooke function>>=
double hooke_handler(double x, void *pdata)
#endif
return k*x;
}
-@
+@
<<hooke structure>>=
typedef struct hooke_param_struct {
double k; /* spring constant in N/m */
} hooke_param_t;
-@
+@
<<hooke tension parameter create/destroy function declarations>>=
extern void *string_create_hooke_param_t(char **param_strings);
extern void destroy_hooke_param_t(void *p);
-@
+@
<<hooke parameter create/destroy functions>>=
hooke_param_t *create_hooke_param_t(double k)
extern int num_hooke_params;
extern char *hooke_param_descriptions[];
extern char hooke_param_string[];
-@
+@
<<hooke tension model globals>>=
int num_hooke_params = 1;
char *hooke_param_descriptions[] = {"spring constant k, N/m"};
char hooke_param_string[]="0.05";
-@
+@
<<hooke tension model getopt>>=
{hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
-@
+@
\subsection{Worm-like chain}
\label{sec.wlc}
<<worm-like chain function>>
<<worm-like chain handler>>
<<worm-like chain create/destroy functions>>
-@
+@
<<worm-like chain tension model declarations>>=
<<worm-like chain handler declaration>>
<<worm-like chain create/destroy function declarations>>
<<worm-like chain tension model global declarations>>
-@
+@
The combination of all unfolded domains is modeled as a worm like chain,
with the force $F_{\text{WLC}}$ approximately given by
$L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
<<worm-like chain function>>=
double wlc(double x, double T, double p, double L)
-{/* N m K m m */
+{/* N m K m m */
double xL=0.0; /* uses global k_B */
assert(x >= 0);
assert(T > 0); assert(p > 0); assert(L > 0);
// k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
}
-@
+@
This model assumes that each unfolded domain has the same persistence length.
<<worm-like chain handler declaration>>=
extern double wlc_handler(double x, void *pdata);
-@
+@
<<worm-like chain handler>>=
double wlc_handler(double x, void *pdata)
#endif
return wlc(x, data->env->T, p, L);
}
-@
+@
<<worm-like chain structure>>=
typedef struct wlc_param_struct {
double p; /* WLC persistence length (m) */
double L; /* WLC contour length (m) */
} wlc_param_t;
-@
+@
<<worm-like chain create/destroy function declarations>>=
extern void *string_create_wlc_param_t(char **param_strings);
extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
-@
+@
<<worm-like chain create/destroy functions>>=
wlc_param_t *create_wlc_param_t(double p, double L)
extern int num_wlc_params;
extern char *wlc_param_descriptions[];
extern char wlc_param_string[];
-@
+@
<<worm-like chain tension model globals>>=
int num_wlc_params = 2;
char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
char wlc_param_string[]="0.39e-9,28.4e-9";
-@
+@
<<worm-like chain tension model getopt>>=
{wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
-@
+@
Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
\subsection{Freely-jointed chain}
<<freely-jointed chain function>>
<<freely-jointed chain handler>>
<<freely-jointed chain create/destroy functions>>
-@
+@
<<freely-jointed chain tension model declarations>>=
<<freely-jointed chain handler declaration>>
<<freely-jointed chain create/destroy function declarations>>
<<freely-jointed chain tension model global declarations>>
-@
+@
The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
The tension of a chain of $N$ such links is given by
//printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
return k_B*T/l * inv_langevin(x/L);
}
-@
+@
<<freely-jointed chain handler declaration>>=
extern double fjc_handler(double x, void *pdata);
-@
+@
<<freely-jointed chain handler>>=
double fjc_handler(double x, void *pdata)
{
#endif
return fjc(x, data->env->T, l, N);
}
-@
+@
<<freely-jointed chain structure>>=
typedef struct fjc_param_struct {
double l; /* FJC link length (m) */
int N; /* FJC number of links */
} fjc_param_t;
-@
+@
<<freely-jointed chain create/destroy function declarations>>=
extern void *string_create_fjc_param_t(char **param_strings);
extern void destroy_fjc_param_t(void *p);
-@
+@
<<freely-jointed chain create/destroy functions>>=
fjc_param_t *create_fjc_param_t(double l, double N)
if (p)
free(p);
}
-@
+@
<<freely-jointed chain tension model global declarations>>=
extern int num_fjc_params;
extern char *fjc_param_descriptions[];
extern char fjc_param_string[];
-@
+@
<<freely-jointed chain tension model globals>>=
int num_fjc_params = 2;
char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
char fjc_param_string[]="0.5e-9,200";
-@
+@
<<freely-jointed chain tension model getopt>>=
{fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
-@
+@
\subsection{Tension model implementation}
<<tension model internal definitions>>
<<tension model globals>>
<<tension model functions>>
-@
+@
<<tension model includes>>=
#include <assert.h> /* assert() */
#include "global.h"
#include "list.h"
#include "tension_balance.h" /* oneD_*() */
-@
+@
<<tension model internal definitions>>=
<<constant tension structure>>
<<hooke structure>>
<<worm-like chain structure>>
<<freely-jointed chain structure>>
-@
+@
<<tension model globals>>=
<<tension handler array global>>
<<hooke tension model globals>>
<<worm-like chain tension model globals>>
<<freely-jointed chain tension model globals>>
-@
+@
<<tension model functions>>=
<<constant tension functions>>
<<hooke functions>>
<<worm-like chain functions>>
<<freely-jointed chain functions>>
-@
+@
\subsection{Utilities}
<<tension model utility getopt functions>>
<<setup>>
<<tension model utility main>>
-@
+@
<<tension model utility main>>=
int main(int argc, char **argv)
(*model->destructor)(params);
return 0;
}
-@
+@
<<tension model utility includes>>=
#include <stdlib.h>
#include "parse.h"
#include "list.h"
#include "tension_model.h"
-@
+@
<<tension model utility definitions>>=
<<version definition>>
<<tension model getopt definitions>>
<<initialize model definition>>
-@
+@
<<tension model utility getopt functions>>=
<<index tension model>>
<<tension model utility help>>
<<tension model utility get options>>
-@
+@
<<tension model utility help>>=
void help(char *prog_name,
printf("\t\tdefaults: %s\n", tension_models[i].params);
}
}
-@
+@
<<tension model utility get options>>=
void get_options(int argc, char **argv, environment_t *env,
}
return;
}
-@
+@
\section{Unfolding rate models}
folded domains exist in a single `folded' state and make instantaneous,
stochastic transitions over a free energy barrier; and Kramers' model, which
treats the unfolding as a thermalized diffusion process.
-We deal with these options like we dealt with the different tension models: by bundling all model-specific
+We deal with these options like we dealt with the different tension models: by bundling all model-specific
parameters into structures, with model specific functions for determining $k$.
<<k func definition>>=
typedef double k_func_t(double F, environment_t *env, void *params);
-@
+@
Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
The interface is defined in [[k_model.h]], the implementation in [[k_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[k_model_utils.c]].
<<null k declarations>>
<<const k declarations>>
<<bell k declarations>>
+<<kbell k declarations>>
<<kramers k declarations>>
<<kramers integ k declarations>>
#endif /* K_MODEL_H */
-@
+@
<<k model module makefile lines>>=
NW_SAWSIM_MODS += k_model
CHECK_BINS += check_k_model
check_k_model_MODS = parse string_eval
-@
+@
[[check_k_model]] also depends on the parse module.
<<k model utils makefile lines>>=
gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
clean_k_model_utils :
rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
-@
+@
\subsection{Null}
\label{sec.null_k}
For domains that are never folded.
<<null k declarations>>=
-@
+@
<<null k functions>>=
-@
+@
<<null k globals>>=
-@
+@
<<null k model getopt>>=
{"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
-@
+@
\subsection{Constant rate model}
\label{sec.k_const}
<<const k functions>>=
<<const k function>>
<<const k structure create/destroy functions>>
-@
+@
<<const k declarations>>=
<<const k function declaration>>
<<const k structure create/destroy function declarations>>
<<const k global declarations>>
-@
+@
<<const k structure>>=
typedef struct const_k_param_struct {
double knot; /* transition rate k_0 (frac population per s) */
} const_k_param_t;
-@
+@
<<const k function declaration>>=
double const_k(double F, environment_t *env, void *const_k_params);
-@
+@
<<const k function>>=
double const_k(double F, environment_t *env, void *const_k_params)
{ /* Returns the rate constant k in frac pop per s. */
assert(p->knot > 0);
return p->knot;
}
-@
+@
<<const k structure create/destroy function declarations>>=
void *string_create_const_k_param_t(char **param_strings);
void destroy_const_k_param_t(void *p);
-@
+@
<<const k structure create/destroy functions>>=
const_k_param_t *create_const_k_param_t(double knot)
<<const k model getopt>>=
{"const", "constant rate transitions", &const_k, num_const_k_params, const_k_param_descriptions, (char *)const_k_param_string, &string_create_const_k_param_t, &destroy_const_k_param_t}
-@
+@
\subsection{Bell's model}
\label{sec.bell}
-Quantitatively, Bell's model gives $k$ as
+Quantitatively, Bell's model gives $k$ as
$$
k = k_0 \cdot e^{F dx / k_B T} \;,
$$
<<bell k functions>>=
<<bell k function>>
<<bell k structure create/destroy functions>>
-@
+@
<<bell k declarations>>=
<<bell k function declaration>>
<<bell k structure create/destroy function declarations>>
<<bell k global declarations>>
-@
+@
<<bell k structure>>=
typedef struct bell_param_struct {
double knot; /* transition rate k_0 (frac population per s) */
double dx; /* `distance to transition state' (nm) */
} bell_param_t;
-@
+@
<<bell k function declaration>>=
double bell_k(double F, environment_t *env, void *bell_params);
-@
+@
<<bell k function>>=
double bell_k(double F, environment_t *env, void *bell_params)
{ /* Returns the rate constant k in frac pop per s.
*/
return p->knot * exp(F*p->dx / (k_B*env->T) );
}
-@
+@
<<bell k structure create/destroy function declarations>>=
void *string_create_bell_param_t(char **param_strings);
void destroy_bell_param_t(void *p);
-@
+@
<<bell k structure create/destroy functions>>=
bell_param_t *create_bell_param_t(double knot, double dx)
<<bell k model getopt>>=
{"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t}
-@
+@
Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
+
+\subsection{Linker stiffness corrected Bell model}
+\label{sec.kbell}
+
+Walton et. al showed that the Bell model constant force approximation
+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...
+
+We adjust Bell's model to
+$$
+ k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
+$$
+where $k_0$ is the unforced unfolding rate,
+$F$ is the applied unfolding force,
+$\kappa$ is the effective spring constant,
+$dx$ is the distance to the transition state, and
+$k_B T$ is the thermal energy\citep{schlierf06}.
+
+<<kbell k functions>>=
+<<kbell k function>>
+<<kbell k structure create/destroy functions>>
+@
+
+<<kbell k declarations>>=
+<<kbell k function declaration>>
+<<kbell k structure create/destroy function declarations>>
+<<kbell k global declarations>>
+
+@
+
+<<kbell k structure>>=
+typedef struct kbell_param_struct {
+ double knot; /* transition rate k_0 (frac population per s) */
+ double dx; /* `distance to transition state' (nm) */
+} kbell_param_t;
+@
+
+<<kbell k function declaration>>=
+double kbell_k(double F, environment_t *env, void *kbell_params);
+@
+<<kbell k function>>=
+double kbell_k(double F, environment_t *env, void *kbell_params)
+{ /* Returns the rate constant k in frac pop per s.
+ * Takes F in N, T in K, knot in frac pop per s, and dx in m.
+ * uses global k_B in J/K */
+ kbell_param_t *p = (kbell_param_t *) kbell_params;
+ assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
+ assert(p != NULL);
+ assert(p->knot > 0); assert(p->dx >= 0);
+ return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
+}
+@
+
+<<kbell k structure create/destroy function declarations>>=
+void *string_create_kbell_param_t(char **param_strings);
+void destroy_kbell_param_t(void *p);
+@
+
+<<kbell k structure create/destroy functions>>=
+kbell_param_t *create_kbell_param_t(double knot, double dx)
+{
+ kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
+ assert(ret != NULL);
+ ret->knot = knot;
+ ret->dx = dx;
+ return ret;
+}
+
+void *string_create_kbell_param_t(char **param_strings)
+{
+ return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
+}
+
+void destroy_kbell_param_t(void *p /* kbell_param_t * */)
+{
+ if (p)
+ free(p);
+}
+@
+
+<<kbell k global declarations>>=
+extern int num_kbell_params;
+extern char *kbell_param_descriptions[];
+extern char kbell_param_string[];
+@
+
+<<kbell k globals>>=
+int num_kbell_params = 2;
+char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
+char kbell_param_string[]="3.3e-4,0.25e-9";
+@
+
+<<kbell k model getopt>>=
+{"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t}
+@
+Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
+
+
\subsection{Kramer's model}
\label{sec.kramers}
<<kramers k functions>>=
<<kramers k function>>
<<kramers k structure create/destroy functions>>
-@
+@
<<kramers k declarations>>=
<<kramers k function declaration>>
<<kramers k structure create/destroy function declarations>>
<<kramers k global declarations>>
-@
+@
<<kramers k structure>>=
//typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
//double *E_params; /* parametrize the energy landscape */
//int n_E_params; /* length of E_params array */
} kramers_param_t;
-@
+@
<<kramers k function declaration>>=
extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
extern double kramers_k(double F, environment_t *env, void *kramers_params);
-@
+@
<<kramers k function>>=
double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
{
//fprintf(stderr,"D = %g, lc = %g, lts = %g, Eb = %g, kbT = %g, k = %g\n", p->D, lc, lts, Eb, kbT, p->D/(lc*lts) * exp(-Eb/kbT));
return p->D/(lc*lts) * exp(-Eb/kbT);
}
-@
+@
<<kramers k structure create/destroy function declarations>>=
void *string_create_kramers_param_t(char **param_strings);
void destroy_kramers_param_t(void *p);
-@
+@
<<kramers k structure create/destroy functions>>=
kramers_param_t *create_kramers_param_t(double D,
For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
Schlierf and Rief used a cubic-spline interpolation routine and the double integral form of Kramers' theory, so we get to pick an actual function to approximate the energy landscape.
We follow \citet{shillcock98} and use
-$$
+$$
E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
\p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
-F x
<<kramers k model getopt>>=
{"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python. For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'. For energies, the position 'x' is also defined. Functions must evaluate to a float representing their output and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t}
-@
+@
\subsection{Kramer's model (natural cubic splines)}
\label{sec.kramers_integ}
<<kramers integ B k functions>>
<<kramers integ k function>>
<<kramers integ k structure create/destroy functions>>
-@
+@
<<kramers integ k declarations>>=
<<kramers integ k function declaration>>
<<kramers integ k structure create/destroy function declarations>>
<<kramers integ k global declarations>>
-@
+@
<<kramers integ k structure>>=
typedef double func_t(double x, void *params);
double F; /* for passing into GSL evaluations */
environment_t *env;
} kramers_integ_param_t;
-@
+@
<<spline param structure>>=
typedef struct spline_param_struct {
gsl_spline *spline; /* GSL spline parameters */
gsl_interp_accel *acc; /* GSL spline acceleration data */
} spline_param_t;
-@
+@
We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
$$
//fprintf(stderr, "integralA = %g\n", result);
return result;
}
-@
+@
$$
\text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
//fprintf(stderr, "integralB = %g\n", result);
return result;
}
-@
+@
With the integrals taken care of, evaluating $k$ is simply
$$
$$
<<kramers integ k function declaration>>=
double kramers_integ_k(double F, environment_t *env, void *kramers_params);
-@
+@
<<kramers integ k function>>=
double kramers_integ_k(double F, environment_t *env, void *kramers_params)
{ /* Returns the rate constant k in frac pop per s. Takes F in N. */
kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
return p->D/kramers_integ_k_integralB(F, env, kramers_params);
}
-@
+@
<<kramers integ k structure create/destroy function declarations>>=
void *string_create_kramers_integ_param_t(char **param_strings);
void destroy_kramers_integ_param_t(void *p);
-@
+@
<<kramers integ k structure create/destroy functions>>=
kramers_integ_param_t *create_kramers_integ_param_t(double D,
<<spline functions>>=
<<spline function>>
<<create/destroy spline>>
-@
+@
<<spline function>>=
double spline_eval(double x, void *spline_params)
spline_param_t *p = (spline_param_t *)(spline_params);
return gsl_spline_eval(p->spline, x, p->acc);
}
-@
+@
<<create/destroy spline>>=
spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
free(p);
}
}
-@
+@
<<kramers integ k global declarations>>=
extern int num_kramers_integ_params;
<<kramers integ k model getopt>>=
{"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t}
-@
+@
Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
<<k model internal definitions>>
<<k model globals>>
<<k model functions>>
-@
+@
<<k model includes>>=
#include <assert.h> /* assert() */
#include <gsl/gsl_spline.h>
#include "global.h"
#include "parse.h"
-@
+@
<<k model internal definitions>>=
<<const k structure>>
<<bell k structure>>
+<<kbell k structure>>
<<kramers k structure>>
<<kramers integ k structure>>
<<spline param structure>>
-@
+@
<<k model globals>>=
<<null k globals>>
<<const k globals>>
<<bell k globals>>
+<<kbell k globals>>
<<kramers k globals>>
<<kramers integ k globals>>
-@
+@
<<k model functions>>=
<<null k functions>>
<<const k functions>>
<<bell k functions>>
+<<kbell k functions>>
<<kramers k functions>>
<<kramers integ k functions>>
-@
+@
\subsection{Unfolding model unit tests}
<<model globals>>
<<k model test suite>>
<<main check program>>
-@
+@
<<k model check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
<<check includes>>
#include "global.h"
#include "k_model.h"
-@
+@
<<k model test suite>>=
<<const k tests>>
<<bell k tests>>
<<k model suite definition>>
-@
+@
<<k model suite definition>>=
Suite *test_suite (void)
<<bell k test case adds>>
return s;
}
-@
+@
\subsubsection{Constant}
<<const k test case defs>>=
TCase *tc_const_k = tcase_create("Constant k");
-@
+@
<<const k test case adds>>=
tcase_add_test(tc_const_k, test_const_k_create_destroy);
tcase_add_test(tc_const_k, test_const_k_over_range);
suite_add_tcase(s, tc_const_k);
-@
+@
<<const k tests>>=
int i;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_const_k_param_t(params);
+ p = string_create_const_k_param_t(params);
destroy_const_k_param_t(p);
}
}
env.T = T;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_const_k_param_t(params);
+ p = string_create_const_k_param_t(params);
for ( F=Fm; F<FM; F+=dF ) {
fail_unless(const_k(F, &env, p)==atof(knot[i]),
"const_k(%g, %g, &{%s}) = %g != %s",
}
}
END_TEST
-@
+@
\subsubsection{Bell}
<<bell k test case defs>>=
TCase *tc_bell = tcase_create("Bell's k");
-@
+@
<<bell k test case adds>>=
tcase_add_test(tc_bell, test_bell_k_create_destroy);
tcase_add_test(tc_bell, test_bell_k_at_zero);
tcase_add_test(tc_bell, test_bell_k_at_one);
suite_add_tcase(s, tc_bell);
-@
+@
<<bell k tests>>=
START_TEST(test_bell_k_create_destroy)
int i;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_bell_param_t(params);
+ p = string_create_bell_param_t(params);
destroy_bell_param_t(p);
}
}
env.T = T;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_bell_param_t(params);
+ p = string_create_bell_param_t(params);
fail_unless(bell_k(F, &env, p)==atof(knot[i]),
"bell_k(%g, %g, &{%s,%s}) = %g != %s",
F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
env.T = T;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_bell_param_t(params);
+ p = string_create_bell_param_t(params);
CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
//fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
// "bell_k(%g, %g, &{%s,%s}) = %g != %g",
}
}
END_TEST
-@
+@
<<kramers k tests>>=
-@
+@
<<kramers k test case def>>=
-@
+@
<<kramers k test case add>>=
-@
+@
<<k model function tests>>=
<<check relative error>>
fail_unless(F = k*x, NULL);
}
END_TEST
-@
+@
\subsection{Utilities}
<<k model utility getopt functions>>
<<k model utility multi model E>>
<<k model utility main>>
-@
+@
<<k model utility main>>=
int main(int argc, char **argv)
(*model->destructor)(params);
return 0;
}
-@
+@
<<k model utility multi model E>>=
double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
double E;
if (k == kramers_integ_k)
E = (*p->E)(x, p->E_params);
- else
+ else
E = kramers_E(0, env, model_params, x);
return E;
}
-@
+@
<<k model utility includes>>=
#include <stdlib.h>
#include "parse.h"
#include "string_eval.h"
#include "k_model.h"
-@
+@
<<k model utility definitions>>=
<<version definition>>
<<kramers k structure>>
<<kramers integ k structure>>
-@
+@
<<k model utility getopt functions>>=
<<index k model>>
<<k model utility help>>
<<k model utility get options>>
-@
+@
<<k model utility help>>=
void help(char *prog_name,
printf(" ...\n");
printf("In special mode, the output depends on the model.\n");
printf("For models defining an energy function E(x), that function is printed\n");
- printf(" #Position (m)\tE (J)\n");
+ printf(" #Distance (m)\tE (J)\n");
printf(" 0.0012\t0.0034\n");
printf(" ...\n");
printf("-m\tChange output to standard mode\n");
printf("\t\tdefaults: %s\n", k_models[i].params);
}
}
-@
+@
<<k model utility get options>>=
void get_options(int argc, char **argv, environment_t *env,
}
return;
}
-@
+@
\section{\LaTeX\ documentation}
noweave -latex -index -delay $< > $@
$(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
cp -f $< $@
-@
+@
and compiled using
<<latex makefile lines>>=
$(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
[[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
<<latex makefile lines>>=
-semi-clean_latex :
+semi-clean_latex :
rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
$(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
$(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
$(BUILD_DIR)/sawsim.bib
clean_latex : semi-clean_latex
rm -f $(DOC_DIR)/sawsim.pdf
-@
+@
\section{Makefile}
BUILD_DIR = build
BIN_DIR = bin
DOC_DIR = doc
-# Note: Cannot use, for example, `./build', because make eats the `./'
-# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
+# Note: Cannot use, for example, `./build', because make eats the `./'
+# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
# Modules (X.c, X.h) defined in the noweb file
-NW_SAWSIM_MODS =
-CHECK_BINS =
+NW_SAWSIM_MODS =
+CHECK_BINS =
# Modules defined outside the noweb file
FREE_SAWSIM_MODS = interp tavl
mkdir $@
# Copy over the source external to sawsim
-# Note: Cannot use, for example, `./build', because make eats the `./'
-# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
+# Note: Cannot use, for example, `./build', because make eats the `./'
+# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
.SECONDEXPANSION:
$(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
| $(BUILD_DIR)
Makefile : $(SRC_DIR)/sawsim.nw
notangle -Rmakefile $< | sed 's/ /\t/' > $@
-@
+@
This is fairly self-explanatory. We compile a dynamically linked
version ([[sawsim]]) and a statically linked version
([[sawsim_static]]). The static version is about 9 times as big, but
x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
F &= k_1 x_1 = k_\text{eff} x \\
- k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
+ k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
= \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
\end{align}
\end{proof}
projects / sawsim.git / commitdiff