python - Writing array to vrt

I would like to write an array from a single band image to a 3 band vrt where each band has the ColorInterpretation() tag of red, green, blue respectively.

My code below;

from osgeo import gdal
import os

#IM = "/path/to/image.tif"
IM = im

### read image ###

ds = gdal.Open(IM)
X = ds.RasterXSize
Y = ds.RasterYSize
band = ds.GetRasterBand(1).ReadAsArray()

### write to 3 bands ###
driver = gdal.GetDriverByName("vrt")
outPath = os.path.join(os.path.split(IM)[0], "test_image.vrt")
outIM = driver.Create(outPath, X, Y, 0, gdal.GDT_Int16)
for i in range(1, 4):

    outIM.AddBand()
    outIM.GetRasterBand(i).SetRasterColorInterpretation(2 + i)
    outIM.GetRasterBand(i).WriteArray(band)
    print outIM.GetRasterBand(i).ReadAsArray()
    print outIM.GetRasterBand(i).GetRasterColorInterpretation()
outIM = None

The output is:

[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]

 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
3
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 

 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
4
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]

 [0 0 0 ..., 0 0 0]]
5

So setting the ColorInterpretation() (unlike when using GTiff, see here), but all values are zero. Why?!

#### UPDATE ####

I have been running this through ipython and just noticed this error on the command line:

Error 1: Writing through VRTSourcedRasterBand is not supported