R代码小小的改进，请R语言朋友们帮忙

亲爱的语言同仁们:
我是一个新R-使用者。现在向大家请教一个问题。我用R建立了一个关于流行病学传播动力学模型SEIR，但是我现在需要对模型中的一个变量加以“时间控制”，我不知道如何修改我的代码，希望得到您的帮助。谢谢。
SEIR模型是一种常用的传染病传播动力学模型，它把人群按照疾病的进展过程分为四种状态，即易感状态(Susceptible)，暴露状态(Exposed)，感染状态(Infectious) 和恢复状态(Recovered)。易感个体在未感染病原体之前为易感状态(S)，当易感者接触到病原体后就转变为暴露状态(E)，随着疾病的进展，暴露者经过疾病的潜伏期后出现了相关的症状，并能将病原菌排出体外感染其他易感者，这时就进入了感染状态(I)，接下来经过疾病的传染期后，感染者随即转变为恢复状态(R)。以上每个状态都在一个给定的概率下向下一个状态进行转化，即易感者(S)暴露后转化为暴露者(E)的概率为β，暴露者(E)经过潜伏期转化为感染者(I)的概率为ω，而感染者(I)经过传染期后恢复(R)的概率为γ。根据以上的模式，疾病随时间的发展过程可以用以下的微分方程进行建模计算。
传播动力学模型SEIR微分方程：


其中，N为t时刻总人口数，S(t)、E(t)、I(t)、R(t)为t时刻四种健康状态下的人口数或人口比例。ω为暴露者向感染者转化的概率，γ为感染者向康复者转化的概率，理论上二者分别等于疾病平均潜伏期和平均传染期的倒数。β为易感者向暴露者转化的概率，反应了疾病的传播能力。该参数通常采用基本再生数R_0进行表示，基本再生数指的是平均一个病例能够感染的二代病例数。R_0参数值可通过人群调查或查阅文献获得。


R语言代码：
#example of very simple SEIR model#####
library (deSolve)
seir_model = function (current_timepoint, state_values, parameters)
{
  # create state variables (local variables)
  S = state_values [1]        # susceptibles
  E = state_values [2]        # exposed
  I = state_values [3]        # infectious
  R = state_values [4]        # recovered
  
  with (
    as.list (parameters),     # variable names within parameters can be used
         {
           # compute derivatives
           dS = (-beta * S * I)
           dE = (beta * S * I) - (delta * E)
           dI = (delta * E) - (gamma * I)
           dR = (gamma * I)
           
           # combine results
           results = c (dS, dE, dI, dR)
           list (results)
         }
    )
}
contact_rate = 10                   # number of contacts per day
transmission_probability = 0.07       # transmission probability
infectious_period = 5                 # infectious period
latent_period = 21                    # latent period
beta_value = contact_rate * transmission_probability
gamma_value = 1 / infectious_period
delta_value = 1 / latent_period
Ro = beta_value / gamma_value
parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)
V = 5000        # exposed hosts
X = 9999        # susceptible hosts
Y = 1           # infectious hosts
Z = 0           # recovered hosts

N = V + X + Y + Z
initial_values = c (S = X/N, E = V/N, I = Y/N, R = Z/N)
timepoints = seq (0, 50, by=1)
output = lsoda (initial_values, timepoints, seir_model, parameter_list)


因为我的研究是对于学校，学校每周六，周末放假，这个时候衡量疾病传染力大小的变量β就成为了0值，现在我想对变量β 施加一个时间控制变量，通过调控该时间变量的状态，可以体现学校节假日规律性停课的变化状态，其中，k表示学校停课的日期，当监测日期t为节假日学校停课时,那么我们定义疾病的传播率β变为0。


我不知道如何修改我的代码，希望有人可以帮助我改进一下代码。谢谢。