\documentclass[12pt]{article}
\usepackage{fancyhdr,fancybox,amssymb,epsfig,amsmath,latexsym,bm}
\usepackage{enumerate}
\usepackage{tikz}
\usepackage[all]{xy}
\usepackage{stackrel}
\usepackage[]{hyperref}
\usetikzlibrary{arrows}
\usetikzlibrary{decorations.markings}

\tikzset{big arrow/.style={decoration={markings,mark=at position 0.6
      with {\arrow[scale=2]{>}}}, postaction={decorate}, shorten >=1pt}}

% Code to make an augmented matrix (Usage: \begin{bmatrix}[ccc|c])
\makeatletter
\renewcommand*\env@matrix[1][*\c@MaxMatrixCols c]{%
  \hskip -\arraycolsep
  \let\@ifnextchar\new@ifnextchar
  \array{#1}}
\makeatother

\topmargin=-.5in
\headsep=0.2in
\oddsidemargin=0in
\textwidth=6.7in
\textheight=9.2in
\footskip=.5in

% Change this line to \def\bold{\mathbf} for ``upright boldface.''
\def\bold{\bm}

\newcommand{\F}{\mathbb{F}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\def\0{\bold{0}}
\def\A{\bold{A}}
\def\B{\bold{B}}
\def\C{\bold{C}}
\def\I{\bold{I}}
\def\M{\bold{M}}
\def\P{\bold{P}}
\def\Q{\bold{Q}}
\def\a{\bold{a}}
\def\b{\bold{b}}
\def\c{\bold{c}}
\def\e{\bold{e}}
\def\p{\bold{p}}
\def\q{\bold{q}}
\def\t{\bold{t}}
\def\u{\bold{u}}
\def\v{\bold{v}}
\def\x{\bold{x}}
\def\<{\langle}
\def\>{\rangle}
\def\hat{\widehat}

\def\And{\wedge}
\def\Or{\vee}
%\def\Not{\neg}
\newcommand{\Not}[1]{\overline #1}



\DeclareMathOperator{\OR}{OR}
\DeclareMathOperator{\NOR}{NOR}
\DeclareMathOperator{\Span}{Span}
\DeclareMathOperator{\rank}{rank}
\newcommand{\vv}[2]{\begin{bmatrix} #1 \\ #2 \end{bmatrix}}
\newcommand{\vvv}[3]{\begin{bmatrix} #1 \\ #2 \\ #3 \end{bmatrix}}

\newcommand{\ds}{\displaystyle}

\renewcommand{\footrulewidth}{1pt}
\newenvironment{proof}[1][Proof.]{\smallskip\begin{trivlist}
    \item[\hskip \labelsep {\sffamily #1}]}{\qed\end{trivlist}\bigskip}
\newenvironment{sol}[1][Solution.]{\smallskip\begin{trivlist}
    \item[\hskip \labelsep {\sffamily #1}]}{\qed\end{trivlist}\bigskip}
\newenvironment{fminipage}
{\setlength{\fboxsep}{15pt}\begin{Sbox}\begin{minipage}}
{\end{minipage}\end{Sbox}\fbox{\TheSbox}}
\newcommand{\qed}{\hfill \ensuremath{\Box}}
\newcommand{\disp}{\displaystyle}


\pagestyle{fancy} \lhead{{\sf Homework 8 $|$ Due March 18 (Friday)
}} \rhead{\thepage} \cfoot{{\sf Math 4500 $|$ Mathematical Modeling
    $|$ Spring 2025 $|$ M.~Macauley}}


\begin{document}

%$\;$

%\noindent \emph{Read}: Chapter 4.1: Global dynamics emerging from
%local interactions: agent-based modeling for the life sciences. By
%D.~Gammack, Elsa.~Schaefer, and H.~Gaff, pp.~105--113.


\begin{enumerate}
  
\item \textbf{Basic Boolean models}. Consider a Boolean model of the
  \emph{lac} operon, based on five variables: mRNA ($M$),
  $\beta$-galactosidase ($B$), \emph{lac} permease ($P$), allolactose
  ($A$), and intracellular lactose ($L$), and the following transition
  functions:
  \[
  (f_m,f_B,f_P,f_A,f_L)=\big(A,\;M,\;M,\;A\vee(L\wedge B),\;
  P\vee(L\wedge\overline{B})\big)
  \]

  \begin{enumerate}
\item What other assumptions are made in this model? (E.g., presence
  or absense of extracellular sugars, time-steps, degredation, etc.)
    
  \item As we saw in class, the dynamics do not accurately reflect the
    behavior of the biological system it is meant to model. Therefore,
    something is wrong. Four of these functions are reasonable;
    justify each one in a single well-written sentence.
    
  \item Explain why one of these functions does not accurately reflect
    the underlying biology and/or the model assumptions. Propose a
    modification, aimed at eliminating the biologically infeasible
    fixed point, $(0,0,0,1,0)$. Give a rationale for your modification and
    specify the biological mechanism or model assumptions that justify
    the change.
  
  \item Use Macaulay2 to convert these Boolean functions into
    polynomials over $\F_2=\{0,1\}$. Then find the fixed points using a
    Gr\"obner basis of a particular polynomial ideal.

  \item Draw the wiring diagram of the old and new Boolean model.

  \item Use the BoolNet package in RStudio to find the attractors and
    plot the state space. Include a print-out or screenshot.
    
  \item Repeat the previous step, but with Cyclone. 

  \item Justify why the long-term behavior of the system (fixed
    points) agrees with what we should expect biologically.
    
  \end{enumerate}

  %%=================================================================

\item \textbf{Fixed points}. The following is a model of the tryptophanase (\emph{tna}) operon in E. coli:
  \[
  \begin{array}{llll}
    f_A=M\And\Not{\gamma} \qquad &&& f_W=\omega_e\And B \\
    f_B=M &&& f_{W_m}=(\omega_{em}\And B)\Or \omega_e\Or W \\
    f_C=\Not{\gamma} &&& f_\gamma=\gamma \\ 
    f_M=C\And\Not{P} &&& f_{\omega_e}=\omega_e \\
    f_P=\Not{W}\And\Not{W_m} &&& f_{\omega_{em}}=\omega_{em} \\
  \end{array}
  \]
  The parameters are encoded by ``frozen'' variables: $\gamma$
  is glucose, $\omega_e$ is high levels of extracellular tryptophan,
  and $\omega_{em}$ represents (at least) medium levels.
  \begin{enumerate}
  \item Use Cyclone to find the attractors. Code is provided on the
    course website.
  \item Fix $(\gamma,\omega_e,\omega_{em})=(0,0,1)$ as constants to
    get a $7$-variable model, which assumes no glucose and medium
    levels of extracellular tryptophan. Use BoolNet to find the
    attactors under an asynchronous update, and summarize your
    findings. This can be done with the command
\begin{verbatim}
  >  getAttractors(tnaModel,type="asynchronous")
\end{verbatim}
How does this compare to the synchronous state space?
\item Explain how/why this model exhibits bistability.
  \end{enumerate}
  
  %%=================================================================

\item \textbf{Decay and degredation}. Recall the $3$-variable ODE model of the
  \emph{lac} operon proposed by Yildirim and Mackey in 2004, where
  $M(t)=$ mRNA, $B(t)=$ $\beta$-galactosidase, and $A(t)=$ allolactose
  (concentrations), respectively.
  \begin{align*}
    \frac{dM}{dt}&=
    \alpha_M\frac{1+K_1(e^{-\mu\tau_M}A_{\tau_M})^n}
          {K+K_1(e^{-\mu\tau_M}A_{\tau_M})^n}
    -\widetilde{\gamma_M}M
    \\ \frac{dB}{dt}&=\alpha_Be^{-\mu\tau_B}M_{\tau_B}-\widetilde{\gamma_B}B
    \\ \frac{dA}{dt}&=
    \alpha_AB\frac{L}{K_L+L}-\beta_AB\frac{A}{K_A+A}-\widetilde{\gamma_A}A
  \end{align*}
  Suppose the exponential decay constants are estimated from the
  literature to be $\widetilde{\gamma_M}=.441$,
  $\widetilde{\gamma_B}=.031$, and
  $\widetilde{\gamma_A}=.55$.
  %%                                                                           
  \begin{enumerate}
  \item Compute the half life for $M$, $B$, and $A$.
  \item Justify the following Boolean model by explaining the logical
    expression defining each transition function:
    \[
    \begin{array}{ll}
    f_M=A & f_{B^\downarrow}=\overline{M}\And B
    \\ f_B=M\Or\left(B\And\overline{B}_1^\downarrow\right)\hspace{10mm} &
    f_A=(B\And L_m)\Or L
  \end{array}
  \]
  What is a reasonble assumption of the approximate timestep assumed by this model? (There are multiple feasible solutions, but you need to justify your answer.)
  \item Find the fixed points of this model using computational algebra.
    Use the variable order $(M,B,B_1^\downarrow,A)=(x_1,x_2,x_3,x_4)$, and include
    your code from Macaulay2.
  \item Does this model exhibit bistability? Justify your answer. \medskip
\end{enumerate}
  
  %%=================================================================

\item \textbf{Reduction}. Consider the following Boolean network model of the \emph{lac} operon.
  \[
  \def\arraystretch{1.4}
  \begin{array}{llll}
    f_1=x_4\And\Not{x_5}\And\Not{x_6}\qquad &&& f_8=x_9\Or x_{10} \\
    f_2=x_1 \qquad &&& f_9=x_3\And\Not{G_e}\And L_e \\
    f_3=x_1 \qquad &&& f_{10}=\Not{G_e}\And\big(L_e\Or(x_3\And L_{em})\big) \\
    f_4=\Not{G_e} \qquad &&& f_{G_e}=G_e\\
    f_5=\Not{x_7}\And\overline{x_8} &&& f_{L_e}=L_e \\
    f_6=x_5\Or(\Not{x_7}\And\Not{x_8}) &&& f_{L_{em}}=L_{em} \\
    f_7=x_2\And x_9 &&&
  \end{array}
  \]
  Here, the variables represent
  \[
  (x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10})=(M,B,P,C,R,R_m,A,A_m,L,L_m).
  \]
  \begin{enumerate}
  \item Use Macaulay2 to reduce this Boolean network, as much as
  possible.
\item Draw the wiring diagram of the reduced network. Find its fixed
  point(s) and use these to determine the fixed point(s) of the
  original network by back-substitution.
  \end{enumerate}
  
  
\end{enumerate}


\end{document}


