Systems mapping allows the genes and genetic interactions for viral dynamics to be identified by incorporating ODEs into a mapping framework. Consider a segregating population composed of

*n* HIV-infected patients genotyped for a set of molecular markers. These patients were repeated sampled to measure uninfected cells (

*x*), infected cells (

*y*) and viral load (

*v*) in their plasma at a series of time points. If specific genes exist to affect the system (1) in Appendix 1, the parameters that specify the system should be different among genotypes. Genetic mapping uses a mixture model-based likelihood to estimate genotype-specific parameters. This likelihood is expressed as

where x

_{
i
} = (

*x*
_{
i
}(

*t*
_{1}), …,

*x*(

^{
t
}
*T*
_{
i
})) , y

_{
i
} = (

*y*
_{
i
}(

*t*
_{1}), …,

*y*(

^{
t
}
*T*
_{
i
})) and v

_{
i
} = (

*v*
_{
i
}(

*t*
_{1}), …,

*v*
_{
i
}(

^{
t
}
*T*
_{
i
})) are the phenotypic values of

*x*,

*y*, and

*v* for subject

*i* measured at

*T*
_{
i
} time points,

*ω*
_{
j|i
} is the conditional probability of QTL genotype

*j* (

*j* = 1, …,

*J*) given the marker genotype of patient

*i*, and

*f*
_{
j
}(x

_{
i
},y

_{
i
},v

_{
i
}) is a multivariate normal distribution with expected mean vector for patient

*i* that belongs to genotype

*j*,

and covariance matrix for subject

*i,*
with
,
and
being (*T*
_{
i
} × *T*
_{
i
}) covariance matrices of time-dependent *x*, *y* and *v* values, respectively, and elements off-diagonal being a (*T*
_{
i
} × *T*
_{
i
}) systematical covariance matrix between the two variables.

For a natural population, the conditional probability of functional genotype given a marker genotype (

*ω*
_{
j|i
}) is expressed in terms of the linkage disequilibria between different loci [

48]. In systems mapping, we incorporate ODEs (1) of Appendix 1 into mixture model (1) to estimate genotypic means (2) specified by ODE parameters for different genotypes, expressed as (

*λ*
_{
j
}
*d*
_{
j
}
*β*
_{
j
}
*a*
_{
j
}
*k*
_{
j
}
*u*
_{
j
}) for

*j* = 1, …,

*J*. Since

*x*,

*y* and

*v* variables obey dynamic system (1) of Appendix 1, the derivatives of genotypic means can be expressed in a similar way. Let

*g*
_{
kj|i
}(

*t*
*μ*
_{
kj|i
}) denote the genotypic derivative for variable

*k* (

*k* =

*x*,

*y*, or

*z*), i.e.,

We use *μ*
_{
kj|i
} to denote the genotypic mean of variable *j* for individual *i* belonging to genotype *j* at an arbitrary point in a time course. The Runge–Kutta fourth order algorithm can be used to solve the ODEs.

Next, we need to model the covariance structure by using a parsimonious and flexible approach such as an autoregressive, antedependence, autoregressive moving average, or nonparametric and semiparametric approaches. Yap et al. [49] provided a discussion of how to choose a general approach for covariance structure modeling. In likelihood (1), the conditional probabilities of functional genotypes given marker genotypes can be expressed as a function of recombination fractions for an experimental cross population or linkage disequilibria for a natural population [48, 50]. The estimation of the recombination fractions or linkage disequilibria can be implemented with the Expectation-Maximization (EM) algorithm.

To demonstrate the usefulness of systems mapping, we assume a sample of *n* HIV-infected patients drawn from a natural human population at random. The sample is analyzed by systems mapping, leading to the detection of a molecular marker which is associated with a QTL that determines the dynamics of drug resistance in a way described by (2) in Appendix 1. At the QTL detected, there are three genotypes *AA*, *Aa* and *aa*, each with a different set of curve parameters (*λ*, *d*, *β*, *β*
_{
m
}
*, a*, *k*, *k*
_{
m
}, *u, ε*) estimated by systems mapping. We assume that these parameters are estimated as (10, 0.01, 0.005, 0.02, 0.5, 10, 10, 3, 0.0001) for genotype *AA*, (12, 0.01, 0.005, 0.02, 0.6, 8, 8, 3, 0.0001) for genotype *Aa,* and (12, 0.008, 0.005, 0.02, 0.55, 8, 12, 4, 0.0001) for genotype *aa*. Using these estimated values, we draw the curves of drug resistance dynamics for each genotype (Figure 2). Pronounced differences in the form of these curves indicate that the QTL plays an important part in determining the resistance dynamics of drugs used to treat HIV/AIDS.

The model for systems mapping described above can be expanded in two aspects, mathematical and genetic, to better characterize the genetic architecture of viral dynamics. The mathematical expansions are to incorporate the drug resistance model (2), the cell diversity model (3) and the unifying resistance and cell diversity model (4). These expansions allow the functional genes operating at different pathways of viral-host reactions to be identified and mapped, making system mapping more clinically feasible and meaningful. The genetic expansions aim to not only model individual genes from the host or pathogen genome but also characterize epistatic interactions between genes from different genomes. This can be done by expanding the conditional probability of functional genes given marker genotypes *ω*
_{
j|i
} using a framework derived by Li et al. [42].

By formulating and testing novel hypotheses, system mapping can address many basic questions. For example, they are

- 1)
How do DNA variants regulate viral dynamics?

- 2)
How do these genes affect the average life-times of uninfected cells, infected cells, and free virus, respectively?

- 3)
How do genes determine the emergence and progression of drug resistance?

- 4)
Are there specific genes that control the possibility of virus eradication by antiviral drug?

- 5)
How important are gene-gene interactions and genome-genome interactions to the dynamic behavior of viral load with or without treatment?