Hi, I have a typical gene expression matrix along with the samples age and gender. How to create a model matrix to test whether there is any association between age and gene expression is different in male and female samples? 

Is it something like this?

x.r$AGE = as.numeric(x.r$AGE)
x.r$Gender = as.factor(x.r$Gender_name)
design <- model.matrix (~x.r$AGE+x.r$Gender, data=x.r)

# fit the linear model
fit <- lmFit(mydata, design)

For starters, you should tag your post with limma alone, otherwise the limma maintainers don't get notified.

To answer you question; it seems that you are looking to test the interaction between sex and age. This is covered in great detail in Section 9.5 of the limma user's guide. If AGE is real-valued, you will probably have to use a classic interaction model, and then perform a test for whether the interaction term is non-zero.

It would also be a good idea to read Section 9.6 to see how to handle continuous variables like time.